BAB
V R untuk Pengolahan & Analisis Statistik
V.1. Statistika deskriptif V.1.1. Rerata (Mean) Rerata merupakan deskripsi statistika yang menggambarkan tentang nilai rata-rata dari suatu sample. Perhitungan rerata secara matematis adalah sebagai berikut: n
mean = 1 / n * ∑ xi i
Dalam R terdapat fungsi untuk menghitung nilai rerata sampel. Fungsi yang digunakan adalah mean(x) ataupun dengan menggunakan fungsi summary(x). Sebagai contoh, digunakan data dari datapackage yang sudah tersedia di R, pilih salah satu data (misalkan Nile). Kemudian hitung nilai rerata sampel, dengan menuliskan > data() > data(Nile) > Nile Time Series: Start = 1871 End = 1970 Frequency = 1 [1] 1120 1160 [14] 994 1020 [28] 1030 1100 [41] 969 831 [54] 864 862 [67] 897 822 [80] 848 890 [94] 906 901
963 960 774 726 698 1010 744 1170
1210 1180 840 456 845 771 749 912
1160 799 874 824 744 676 838 746
1160 958 694 702 796 649 1050 919
813 1140 940 1120 1040 846 918 718
1230 1100 833 1100 759 812 986 714
1370 1210 701 832 781 742 797 740
1140 1150 916 764 865 801 923
995 1250 692 821 845 1040 975
935 1260 1020 768 944 860 815
1110 1220 1050 845 984 874 1020
47
> mean(Nile) [1] 919.35 > summary(Nile) Min. 1st Qu. 456.0 798.5
Median 893.5
Mean 919.4
3rd Qu. 1033.0
Max. 1370.0
Hasil tersebut di atas adalah menunjukkan tentang data Nile yang terdiri dari 100 data dengan nilai rerata 919.35. Selain rerata ada pula nilai statistik lainnya, yaitu minimum, kuartil pertama, nilai tengah (median), kuartil ke tiga dan maksimum. Nilai-nilai tersebut ditampilkan dengan menjalankan fungsi summary(). V.1.2. Nilai Tengah (median)
Seperti halnya dengan rerata, R juga menyediakan fungsi untuk mencari nilai tengah (median) sampel dengan menuliskan fungsi median(). Sebagai contoh, dapat digunakan data seperti mencari nilai rerata sebelumnya untuk menghitung nilai tengah (median): > median(Nile) [1] 893.5
V.1.3. Nilai minimum dan maksimum
R menyediakan fasilitas untuk mencari nilai minimum dan maksimum suatu data, yaitu dengan digunakan perintah min() dan max() > min(x) > max(x)
# untuk nilai minimum, dan # untuk nilai maksimum
V. 2. Grafik
R dilengkapi dengan fasilitas untuk visualisasi statistik dalam bentuk grafik, baik statistik, kontur, map, dll. Sistem grafik di R terdiri dari dua sistem: sistem (dasar/default) yang terdapat dalam paket graphics dan sistem trellis yang terdapat dalam paket lattice. Grafik di R dapat diatur sesuai keperluan. Untuk melihat lebih detil mengenai fitur grafik dalam R, pilih menu Help di menu R kemudian pilih Manual (dalam format pdf) atau HTML help. Untuk mendapatkan gambaran langsung tentang grafik dalam R, dapat dilihat dalam fungsi demo(), dengan menuliskan > > > > >
demo() demo(graphics) demo(image) demo(persp) demo(lattice)
# # # # #
untuk melihat jenis-jenis demo atau atau atau sebelumnya diperlukan load package dan pilih lattice pada menu di windows R anda; atau > demo(package = .packages(all.available = TRUE)) # untuk melihat semua jenis demo yang tersedia
Ket: untuk R versi 2.1.1, penulisan Return dituliskan setelah penulisan fungsi demo()
48
Grafik Dasar (Base)
Pengantar tentang prosedur grafik dapat dilihat di dokumen “Introduction to R” pada menu Help ->Manual (dalam format pdf). Berikut akan disajikan contoh pembuatan plot secara bertahap diawali dengan model standar hingga pengaturan sesuai dengan yang diinginkan (customize). Contoh berikut adalah pembuatan scatterplot untuk petal.length yang dibandingkan dengan petal.width dari dataset iris. Default scatterplot dari dua variabel dihasilkan oleh metod plot.default, yang secara otomatis digunakan oleh perintah plot generik dimana argumennya merupakan dua vektor dengan panjang yang sama seperti berikut ini: > data(iris) > str(iris) ‘data.frame’: 150 obs. of 5 variables: $ Sepal.Length: num 5.1 4.9 4.7 4.6 5 5.4 4.6 5 4.4 4.9 … $ Sepal.Width : num 3.5 3 3.2 3.1 3.6 3.9 3.4 3.4 2.9 3.1 … $ Petal.Length: num 1.4 1.4 1.3 1.5 1.4 1.7 1.4 1.5 1.4 1.5 … $ Petal.Width : num 0.2 0.2 0.2 0.2 0.2 0.4 0.3 0.2 0.2 0.1 … $ Species : Factor w/ 3 levels "setosa","versicolor",..: 1 1 1 1 1 1
Data tersebut di atas menyatakan bahwa data iris terdiri dari 5 variabel dimana setiap variable terdiri dari 150 data observasi. Lima variable tersebut adalah: Sepal.Length, Sepal.Width, Petal.Length, Petal.Width dan Species. > attach(iris) > plot(Petal.Length, Petal.Width)
Hasil dari perintah tersebut terlihat pada Grafik 5.1 (merupakan bentuk grafik secara default). Dalam grafik tersebut, sumbu x dan sumbu y berasal dari argumen pertama (Petal.Length) dan argumen kedua (Petal.Width), dari pernyataan plot(argument_1, argument_2).
49
Gambar 5.1: Scatter plot data variabel Petal Grafik 5.1 di atas dapat dilengkapi untuk menunjukan dependensi argumen dimana sumbu y sebagai variabel dependen. Hal tersebut dilakukan dengan menuliskan perintah seperti berikut, dimana variable terikat terletak sebelah kiri: > plot(Petal.Width ~ Petal.Length,)
Pada Gambar 5.1 bentuk grafik sangat standart, sehingga perlu dilengkapi dengan beberapa keterangan tambahan untuk memperjelas dan mempermudah dalam melakukan interpretasi grafik. Hal ini dapat dilakukan dengan menambahkan fitur warna atau simbol dalam tampilan grafik. Untuk hal tersebut, R mempunyai fasilitas pewarnaan (yaitu dengan argumen col), simbol (dengan argumen pch), ukuran (dengan argumen cex), label/nama sumbu kordinat (dengan argumen xlab dan ylab), judul grafik (dengan argumen main). Beberapa jenis warna yang disediakan dalam R dapat diketahui dengan menggunakan perintah colours(), dimana akan ditampilan daftar warna-warna tersebut. > colours() [1] "white" [4] "antiquewhite1" … [655] "yellow3"
"aliceblue" "antiquewhite2
"antiquewhite" " antiquewhite3"
"yellow4"
"yellowgreen"
sedangkan fungsi palette() dapat digunakan untuk menampilkan daftar warna dalam bilangan numeric: > palette() [1] "black" "red" "green3" "blue" "cyan" "magenta" "yellow" [8] "gray"
50
Selain warna, simbol juga dapat digunakan untuk menampilkan tanda plot. Spesifikasi simbol dapat dilakukan dengan menentukan karakter yang akan digunakan (misalkan asterik ”*”) atau kode integer dari simbol tersebut. Gambar 5.2 menunjukkan simbol dan kodenya. Simbol dengan kode 21-26 memiliki fill (warna latar) yang dispesifikasikan pada argumen bg, spesifikasi warna utama dengan argumen col menspesifikasikan garis border.
Gambar 5.2: Simbol dan kode dalam R Berikut ini merupakan contoh yang menampilkan grafik plot yang dilengkapi dengan pewarnaan, modifikasi ukuran dan pemilihan simbol plotting. Penulisannya adalah sebagai berikut: > > > > > > >
plot(1, 1, xlim=c(1, 7.5), ylim=c(0,5), type="n") points(1:7, rep(4.5, 7), cex=1:7, col=1:7, pch=0:6) text(1:7,rep(3.5, 7), labels=paste(0:6), cex=1:7, col=1:7) points(1:7,rep(2,7), pch=(0:6)+7) # Plot simbol 7 hingga 13 text((1:7)+0.25, rep(2,7), paste((0:6)+7)) # Label dengan bilangan simbol points(1:7,rep(1,7), pch=(0:6)+14) # Plot symbols 14 hingga 20 text((1:7)+0.25, rep(1,7), paste((0:6)+14)) # Labels dengan bilangan simbol
dan akan menghasilkan Gambar 5.3 sebagai berikut:
51
Gambar 5.3: Simbol, Kode dan Warna dalam R Setelah kita mengetahui bentuk dan kode simbol, maka berikut ini akan dilakukan penggunaan simbol, warna dan modifkasi lain untuk menggambarkan plot/grafik dari contoh sebelumnya dengan menuliskan perintah seperti berikut: > + + + +
plot(Petal.Length, Petal.Width, pch=20, cex=1.2, xlab=" Petal length (cm)", ylab="Petal width (cm)", main="Data Anderson Iris", col=c("slateblue","firebrick",darkolivegreen") [as.numeric(Species)])
Data menunjukkan dengan jelas bahwa ukuran spesies berbeda (Sentosa paling kecil, Versicolor menengah, Virginica terbesar) tetapi rasio petal length dan weight sama untuk ketiga ukuran tersebut.
Gambar 5.4: Grafik scatterplot data Iris
52
V.2.1. Histogram
Selain plot, bentuk representasi grafis lainnya yang paling mudah digunakan untuk menggambarkan sebaran data adalah histogram. R menyediakan fasilitas fungsi histogram yang digunakan untuk mengetahui sebaran sampel suatu data. Sebagai catatan: histogram ataupun boxplot, digunakan untuk satu variable. Sebelum kita mencoba untuk menggunakan fasilitas histogram, maka perlu sedikit penjelasan yang berkaitan dengan histogram, yaitu: • histogram digunakan untuk mengestimasi fungsi distribusi probabilitas densitas (probability density function);
f ( x) = lim Pr ob( x − δ < X ≤ x) / δ δ →0
• • • • •
histogram ditentukan pula oleh bin/lebar batang; sumbu- y dalam histogram dapat berupa frekuensi kemunculan atau proporsi; tidak ada estimasi statistik yang dapat dibaca langsung dari histogram, namun dengan histogram kita dapat menduga kemiringan, sifat/behavior pada tail atau ujung kurva, dan outlier data; histogram dapat dibandingkan sebagai suatu distribusi analitik standar.
Selain histogram, R juga menyediakan plot yang fungsinya menyerupai histogram yakni stem-and-leaf plot dengan penulisan: > stem(variabel)
Sebagai contoh, gunakan variabel eruptions dari tabel data faithful, dengan menuliskan: > attach(faithful) > summary(eruptions) Min. 1.600
1st Qu. Median 2.163 4.000
Mean 3.488
3rd Qu. 4.454
Max. 5.100
> fivenum(eruptions) [1] 1.6000 2.1585 4.0000 4.4585 5.1000 > stem(eruptions) 16 18 20 22 24 26 28 30 32 34 36 38
| | | | | | | | | | | |
The decimal point is 1 digit(s) to the left of the | 070355555588 000022233333335577777777888822335777888 00002223378800035778 0002335578023578 00228 23 080 7 2337 250077 0000823577 2333335582225577
53
40 42 44 46 48 50
| | | | | |
0000003357788888002233555577778 03335555778800233333555577778 02222335557780000000023333357778888 0000233357700000023578 00000022335800333 0370
Kita juga dapat melihat sebaran data dalam plot histogram yaitu dengan menggunakan fungsi hist() > hist(eruptions)
Gambar 5.5: Grafik histogram data eruptions Pada Gambar 5.5 di atas, fungsi hist() menggunakan jarak antar batang (disebut bin) cukup besar. Untuk membuat bin lebih kecil, diperlukan tambahan atribut dengan menuliskan: > hist(eruptions, seq(1.6, 5.2, 0.2), prob=TRUE)
Pada perintah di atas, argumen seq(1.6, 5.2, 0.2) adalah histogram menggunakan range dari 1.6 hingga 5.2 dengan lebar bin 0.2. Sehingga tampilan grafik histogram adalah sebagai berikut:
54
Gambar 5.6: Grafik histogram data eruption berdasarkan lebar bin 0.2 Gambar 5.6 di atas menunjukkan lebar batang histogram yang lebih kecil dibanding Gambar 5.5. Apabila ingin ditambahkan garis pada data densitas, maka dapat menggunakan fungsi lines()seperti berikut: > lines(density(eruptions, bw = 0.1))
Keterangan: bw adalah bandwidth (lebar pita), dengan nilainya berdasarkan trial and error.
Gambar 5.7: Grafik histogram data eruptions dengan bw = 0.1 Untuk menampilkan point data aktual digunakan fungsi rug() sbb: > rug(eruptions)
55
Gambar 5.8: Grafik histogram data eruption dengan data aktual Selain histogram, R dapat memplot fungsi distribusi kumulatif empiris dengan menggunakan fungsi ecdf(): > plot(ecdf(eruptions), do.points=FALSE, verticals=TRUE)
Gambar 5.9: Plot 1 ecdf data eruption Distribusi ecdf di atas masih jauh dari standar distribusi yang ada. Untuk itu dapat dicobakan dengan mencocokkan distribusi normal dan “menutupi” fungsi distribusi kumulatif (ecdf) sebelumnya. Penulisannya adalah sebagai berikut: > long <- eruptions[eruptions > 3] > plot(ecdf(long), do.points=FALSE, verticals=TRUE)
56
Gambar 5.10: Plot 2 ecdf data eruption > x <- seq(3, 5.4, 0.01) > lines(x, pnorm(x, mean=mean(long), sd=sqrt(var(long))), lty=3)
Gambar 5.11: Plot 3 ecdf data eruption V.2.2. Q-Q (Quantile – Quantile)
Selain histogram sebagai alat untuk memplot sebaran data suatu variabel adalah Quantile – Quantile (Q-Q) plot. Q-Q plot dapat digunakan untuk memplot variable secara lebih teliti berdasarkan nilai quantile data. Q-Q plot adalah suatu scatter plot yang membandingkan distribusi empiris dengan fitted distribution dalam kaitannya dengan nilai dimensi suatu variabel (misalkan: nilai quantile empiris). Q-Q plot dapat memplot dengan baik jika dataset diperoleh dari populasi yang sudah diketahui. 57
Q-Q plot dalam R dibagi menjadi dua, yaitu: • qqnorm(variabel); untuk menguji goodness of fit dari distribusi Gaussian. qqnnorm() disebut juga sebagai plot probabilitas normal. •
qqplot(variabel); untuk sebarang jenis distribusi
Melanjutkan contoh sebelumnya, untuk data eruptions, penggunaan Q-Q plot dituliskan sbb: > par(pty="s") # mengatur pembuatan daerah gambar berbentuk kotak > qqnorm(long)
Gambar 5.12: Plot1 Q-Q Normal > qqline(long)
Gambar 5.13: Plot 2 Q-Q Normal
58
V.2.3. Boxplot
Selain dua alat untuk menggambarkan grafik untuk satu variable yang sudah dijelaskan sebelumnya, terdapat fasilitas boxplot yang digunakan untuk melihat sebaran data. Berikut adalah penjelasan tentang fitur dasar boxplot: • Berguna untuk membanding banyak kelompok/grup. • Dasarnya menggunakan 3 jenis summary: 3 quartil. • Mudah dalam menampilkan nilai rerata (mean). • Dapat diperluas untuk menampilkan persentil lainnya, terutama pada ujung(tails) suatu distribusi. R menyediakan fitur untuk menampilkan boxplot, dengan menuliskan fungsi boxplot(variabel). Untuk menjelaskan penggunaan fungsi boxplot(), berikut adalah contoh menggambar grafik dengan menggunakan data faithful seperti pada contoh sebelumnya. > boxplot(eruptions)
Perintah di atas akan mendapatkan gambar boxplot dari variabel eruption seperti berikut:
Gambar 5.14: Boxplot data eruption Apabila ingin menambahkan judul gambar, maka penulisannya adalah: > boxplot(eruptions, main=”Plot dengan Boxplot”)
59
Gambar 5.15: Boxplot data eruption dengan nama titel R menyediakan beragam bentuk penyajian grafik plot. Berikut adalah daftar plot grafik dasar yang ada dalam R (beberapa ada yang termasuk dalam instalasi awal dan adapula yang masuk dalam paket lain yang harus didownload dari CRAN): Tabel 5.1: Jenis plot dalam R Jenis Fungsi Plot Keterangan assocplot Plot Asosiasi (Association) barplot Plot Batang (Bar) boxplot Plot Kotak (Box) contour Plot Contour coplot Plot Conditioning dotchart Plot Cleveland Dot filled.contour Plot Level (Contour) fourfoldplot Plot Fourfold hist Histogam image Menampilkan suatu Warna Image matplot Plot Kolom suatu Matriks mosaicplot Plot Mosaic pairs Matriks Scatterplot persp Plot Perspektif plot Plot X-Y Umum stars Plot Star (Spider/Radar) stem Plot Stem-and-Leaf stripchart Plot Scatter 1-D sunflowerplot Plot Scatter Sunflower
60
Gambar 5.16 terdiri dari beberapa contoh tampilan grafik plot, yaitu boxplot, a conditioning plot, pairwise scatterplot, dan star plot, yang kesemuanya mengaplikasikan dataset Anderson iris. Perintah di R untuk menggambarkan grafik-grafik tersebut adalah : > + + + > + > + > +
boxplot(Petal.Length ~ Species, horizontal=T, col="lightblue", boxwex=.5, xlab="Petal length (cm)", ylab="Species", main="BoxPlot Group") coplot(Petal.Width ~ Petal.Length | Species, col=as.numeric(Species), pch=as.numeric(Species)) pairs(iris[,1:4], col=as.numeric(Species), main="Pairwise scatterplot") stars(iris[,1:4], key.loc=c(2,35), mar=c(2, 2, 10, 2), main="Star plot of individuals", frame=T)
Conditioning plot
boxplot
Pairwise scatterplot
Star plot
Gambar 5.16: Contoh grafik plot
61
V.2.4 Grafik Trellis
Sistem grafik trellis dalam R tersedia dalam paket lattice. Model grafik ini khususnya digunakan untuk visualisasi multivariate apabila relasi antara variable berubah bersama beberapa group factor yang disebut sebagai kondisi (conditioning) suatu grafik terhadap factor. Metoda ini menggunakan formula yang similar dengan formula statistic untuk menspesifikasikan variable yang akan diplot serta hubungannya dalam plot. Untuk memudahkan dalam penggunaannya, jenis plot dibagi berdasarkan banyaknya variabel: satu variable (univariate), dua variable (bivariate), tiga variable (trivariate) dan banyak variable (hypervariate). Satu varibel/Univariate
Sebagai salah contoh univariate adalah membuat grafik plot densitas pada keseluruhan data. Pada contoh ini digunakan dataset iris seperti contoh sebelumnya. Berikut adalah metode yang digunakan untuk menampilkan plot densitas: > densityplot(~ Petal.Length, data=iris,main=”Seluruh Spesies”)
Operator ~ tidak memiliki operan di sebelah kiri, karena tidak ada variable terikat (dependent) dalam plot; ini menunjukkan sifat univariate. Petal.Length adalah variable bebas (independent), dan diperoleh plotnya. Gambar 5.17(a) sebelah kiri menunjukkan plot densitas univariate. Pengkondisian dilakukan dengan menambahkan operator “|”, yang dapat dibaca sebagai “pengkondisian pada” (conditioned on) satu (beberapa) variable pada sisi kanan operator, seperti pada berikut ini: > densityplot(~ Petal.Length | Species, data=iris)
Perintah tersebut akan menampilkan satu panel per spesies; seperti yang ditunjukkan pada Gambar 5.17(b). Pada gambar tersebut tampak jelas bahwa distribusi multimodal keseluruhan data set berdasarkan perbedaan distribusi untuk setiap spesies.
(a) (b) Gambar 5.17 Plot densitas triller (a) tanpa dan (b) dengan pengkondisian faktor 62
Jenis plot untuk satu variable (univariate) adalah seperti pada table 5.2 berikut: Tabel 5.2: Jenis plot untuk satu variable Jenis Keterangan assocplot Plot hubungan (association) barchart Plot batang (bar) bwplot Plot box dan whisker densityplot Plot kepadatan Kernel dotplot Plot dot histogram Histogram qqmath Quantile plot distribusi matematis stripplot Scatterplot 1 dimensi
Dua variable (Bivariate)
Salah satu metode membuat plot dua variable (bivariate) adalah xyplot, dimana sumbu y adalah variable terikat dan sumbu x adalah variable bebas; variable tersebut juga dapat dikondisikan terhadap satu atau lebih kelompok faktor seperti pada perintah berikut. > + > +
xyplot(Petal.Width ~ Petal.Length, data=iris, groups=Species, auto.key=T,main=”Seluruh Spesies”) xyplot(Petal.Width ~ Petal.Length | Species, data=iris, groups=Species,main=”Dipisah berdasarkan spesies”)
Gambar 5.18: Scatter plot triller (a) tanpa dan (b) dengan pengkondisian faktor Gambar 5.18 menunjukkan penggunaan argumen group untuk menspesifikasikan perbedaan cara menampilkan grafik (dalam hal ini warna) untuk setiap spesies, dan argument auto.key untuk mendapatkan kunci sederhana terhadap warna yang digunakan untuk setiap spesies. Motode lain yang ada di R untuk membentuk grafik dengan bivariate adalah pada table 5.3 berikut:
63
Tabel 5.3: Jenis plot untuk dua variable Jenis Keterangan qq Plot untuk membandingkan dua distribusi xyplot Plot scatter Tiga variable(Trivariate)
Plot yang paling banyak digunakan untuk trivariate adalah levelplot dan contourplot untuk melakukan plot 2D dari satu variable respon pada dua variable terikat kontinu (misalkan, elevation vs. dua koordinat), metode wireframe untuk suatu versi grafik 3D, dan metode cloud (awan) untuk scatter plot 3D dari tiga variable. Semua dapat dikondisikan pada suatu factor tertentu. Gambar 5.19 menunjukkan contoh yang dihasilkan dari kode berikut: > + + + > > + + + + > + + > + + + + + > > > > >
pl1 <- cloud(Sepal.Length ~ Petal.Length * Petal.Width, groups=Species, data=iris, pch=20, main="Anderson Iris data, all species", screen=list(z=30, x=-60)) data(volcano) pl2 <- wireframe(volcano, shade = TRUE, aspect = c(61/87, 0.4), light.source = c(10, 0, 10), zoom=1.1, box=F, scales=list(draw=F), xlab="", ylab="", zlab="", main="Wireframe plot, Maunga Whau Volcano, Auckland") pl3 <- levelplot(volcano, col.regions=gray(0:16/16), main="Levelplot, Maunga Whau Volcano, Auckland") pl4 <- contourplot(volcano, at=seq(floor(min(volcano)/10)*10, ceiling(max(volcano)/10)*10, by=10), main="Contourplot, Maunga Whau Volcano, Auckland", sub="contour interval 10 m", region=T, col.regions=terrain.colors(100)) print(pl1, split=c(1,1,2,2), more=T) print(pl2, split=c(2,1,2,2), more=T) print(pl3, split=c(1,2,2,2), more=T) print(pl4, split=c(2,2,2,2), more=F) rm(pl1, pl2, pl3, pl4)
Gambar 5.19 menunjukkan hasil dari perintah di atas. Sebagai catatan, data set volcano merupakan matriks elavasi: > str(volcano) num [1:87, 1:61] 100 101 102 103 104 105 105 106 107 108 ...
Metode levelplot menkonversi ke variable respon (nilai z) dan dua predictor, yaitu baris dan kolom matriks (nilai x dan y). Contoh tersebut menunjukkan metode lattice tingkat tinggi yang melakukan pembentukan grafik sendiri. Hasil dari metode levelplot digambar dengan metode print. Metode plot tiga variable adalah seperti Table 5.4:
64
Gambar 5.19: Grafik plot Trellis trivariate Tabel 5.4: Jenis plot untuk tiga variable Jenis Keterangan levelplot Plot level contourplot Plot contour cloud Plot scatter 3 dimensi wireframe Permukaan 3 dimensi (similar dengan plot persp di R)
Lebih dari tiga variabel/Hypervariate
Selain motode plot untuk satu, dua, tiga variable, R juga menyediakan plot grafik untuk lebih dari tiga variable (hypervariate) seperti pada Table 5.5 berikut: Tabel 5.5: Jenis plot untuk lebih dari tiga variable Jenis Keterangan splom Matriks plot scatter parallel Plot koordinat paralel
V.3. Fungsi Distribusi
Fungsi distribusi merupakan salah satu bahasan penting dalam statistika, terutama dalam analisis data. Fungsi distribusi merupakan salah satu alat pendekatan distribusi suatu data. Fungsi distribusi juga berperan dalam menentukan densitas suatu fungsi
65
data. Dalam bab ini akan dibahas fungsi distribusi dan fungsi yang berkaitan dengannya. V.3.1. Jenis fungsi distribusi dalam R
Software R mempunyai koleksi fungsi distribusi standar yang lengkap, yang tersedia dalam paket program R dan dapat ditambah dengan mendownload dalam bentuk paket dari situs R. Fungsi distribusi di R disediakan untuk memfasilitasi fungsi distribusi kumulatif (Cummulative Distributive Function (CDF)) P(X ≤ x),fungsi probabilitas densitas (Probability Density Function (PDF)), dan fungsi kuantil (diberikan q, x lebih kecil sedemikian hingga P(X ≤ x) > q). Berikut adalah tabel distribusi di R.
Tabel 5.6: Jenis fungsi distribusi dalam R Nama Distribusi Nama fungsi di R Argument tambahan beta Beta shape1, shape2, ncp binom Binomial size, prob cauchy Cauchy location, scale chisq Chi-squared df, ncp exp Exponential rate f F df1, df1, ncp gamma Gamma shape, scale geom Geometric Prob hyper Hypergeometric m, n, k lnorm Log-normal meanlog, sdlog logis Logistic location, scale nbinom Binomial negative size, prob norm Normal mean, sd pois Poisson Lambda t t - Student’s df, ncp unif Uniform min, max weibull Weibull shape, scale wilcox Wilcoxon m, n
V.3.2. Fungsi Probabilitas Densitas (Probability Density Function (PDF))
Fungsi probabilitas densitas merupakan salah satu parameter statistic yang digunakan untuk mengetahui probabilitas terhadap suatu factor yang ada dalam sample. Buku ini akan mengawali pembahasan dengan contoh kasus untuk mempermudah pemahaman tentang fungsi densitas seperti berikut ini: Terdapat 16 mahasiswa dipilih secara acak dari populasi dimana 30% adalah wanita. Berapa probabilitas sebanyak nol, satu, dua, …, enam belas dari mahasiswa tersebut yang dipilih adalah wanita?. Untuk menghitung probabilitas tersebut akan digunakan beberapa langkah dalam R seperti dibawah ini.
66
> round(dbinom(0:16, 16, 0.3), 3)
# dbinom artinya d:densitas dan binom:binomial
[1] 0.003 0.023 0.073 0.146 0.204 0.210 0.165 0.101 0.049 0.019 0.006 0.001 [13] 0.000 0.000 0.000 0.000 0.000
Pertama adalah menghitung nilai probabilitas jumlah wanita yang terpilih dari populasi yang dicari tersebut berdasarkan distribusi binomial dengan menggunakan fungsi dbinom(). Nilai probabilitas dari masing-masing kejadian tersebut adalah: Jml Wanita terpilih 0 1 2 3 4 5 6 7 8 9 10 11 12 - 16
Probabilitas 0.003 0.023 0.073 0.146 0.204 0.210 0.165 0.101 0.049 0.019 0.006 0.001 0.000
Total kumulatif dari masing-masing nilai probabilitas tersebut adalah 1 (hal ini menunjukkan bahwa nilai probabilitas tersebut merupakan pdf). Kemudian untuk menggambarkan scatter plot nilai probabilitas tersebut digunakan perintah plot seperti berikut: > plot(dbinom(0:16, 16, 0.3), type =”h”, xlab=” + ylab=”Probabilitas”)
Banyaknya wanita”,
Perintah tersebut menggambarkan grafik plot (Gambar 5.20) dimana sample berdistribusi binomial dimana menghitung probabilitas sukses 0 hingga 16 (0:16) pada suatu populasi yang terdiri dari 16 mahasiswa dimana terdiri dari 0.3 (30%) mahasiswa wanita. Misalkan terdapat kondisi bahwa hanya 2 dari 16 yang terpilih adalah wanita. Berapa probabilitas dua atau kurang dari jumlah wanita dapat terpilih kembali? Penyelesaian masalah tersebut dapat diselesaikan dengan menggunakan fungsi pbiom() seperti perintah program dibawah, dimana argument pertama = 2 yang menyatakan jumlah wanita terpilih, argument kedua adalah 16 yang menyatakan jumlah keseluruhan, argument ketiga menyatakan probabilitas jumlah wanita. Perintah dalam R adalah sebagai berikut: > pbinom(2, 16, 0.3, lower.tail = T) # pbinom, p: fungsi distribusi atau probabilitas, binom: binomial [1] 0.09935968
67
Gambar 5.20: Plot binomial data mahasiswa (Interpretasi: dalam sample random dari 16 orang populasi dengan 30% adalah wanita, terdapat sekitar 10% bahwa dari sample akan terpilih sebanyak nol, satu atau dua adalah wanita. Maka jika kita lihat, 16 orang sample dengan dua atau kurang wanita, maka kita menduga bahwa terjadi diskriminasi terhadap hal tersebut, maka 10% kesempatan yang kita duga tersebut adalah tanpa dasar yang kuat). Contoh berikut menyatakan bahwa dua random sample yang saling bebas yang berdistribusi normal seharusnya tidak berkorelasi. Kita dapat mensimulasikan hal ini berulangkali untuk menaksi koefisien korelasi dimana probabilitas error Type I (yaitu menolak hipotesis yang menyatakan tidak terdapat suatu korelasi) adalah 10%. Pertama, kita akan menuliskan: > size <- 20 > sample.1 <- rnorm(size) # rnorm; r: deviasi random; norm: distribusi normal > sample.2 <- rnorm(size) > cor.test(sample.1, sample.2) Pearson’s product-moment correlation data: sample.1 and sample.2 t = 0.4394, df = 18, p-value = 0.6656 alternative hypothesis: true correlation is not equal to 0 95 percent confidence interval: -0.3557135 0.5217580 sample estimates: cor 0.1030244 > > > >
plot(sample.1, sample.2) abline(h=0, lty=2) abline(v=0, lty=2) title(paste(“r=”, round(cor(sample.1, sample.2), 3)))
68
Gambar 5.21: Plot variabel sample1 dan sample 2 Gambar 5.21 menunjukkan bahwa 2 variabel random sample.1 dan sample.2 adalah saling bebas. Sebagai catatan, karena sifat kerandomannya, hasil yang akan anda peroleh dari hasil percobaan yang dilakukan tentu akan berbeda dengan apa yang disajikan dalam buku ini. Kemudian, untuk perhitungan yang lebih besar lagi, yakni mencari nilai-nilai statistiknya, perlu mendefinisikan fungsi secara tersendiri seperti berikut ini: > cor.2 <- function(size) { + sample.1 <- rnorm(size) + sample.2 <- rnorm(size) + cor(sample.1, sample.2) +} > results <- NULL > length <- 1000 > for (i in 1:length) { + results[i] <- cor.2(20) +} > hist(results, xlab= “r”, main= “Koefisien Korelasi”) > (paste(“5% nilai r lebih negative dari”, round(sort(results)[length/20], 3))) [1] "5% nilai r lebih negative dari -0.386" > (paste(“5% nilai r lebih positif dari”, round(sort(results)[length - length/20], 3))) [1] "5% nilai r lebih positif dari 0.35"
69
Gambar 5.22: Histogram koefisien korelasi Sepertinya, terlihat suatu interval antara -0, 36 …. +0,37 yang melingkupi 90% dari koefisien korelasi r sample untuk ukuran sample 20. Sekali lagi hasilnya mungkin sedikit berbeda antara yang anda lakukan dengan yang dikerjakan dalam buku ini. Suatu plot perspektif untuk PDF - distribusi normal bivariat: > > > > > + +
range <- seq (-3, 3, by= 0.1) random.1 <- dnorm(range, 0, 1) random.2 <- dnorm(range, 0, 0.7) mesh <- outer(random.1, random.2, function(x, y) x * y) persp(range, range, mesh, phi = 20, theta = 30, expand = 0.8, xlab = “N(0, 1)”, ylab = “N(0, 0.7)”,zlab = “PDF”, main = ”Contoh PDF untuk Distribusi Normal Bivariat”)
Gambar 5.23: Grafik 3D untuk pdf normal bivariat
70
Baris perintah di atas adalah salah satu contoh pembahasan mengenai fungsi probabilitas densitas (probability density function (PDF)) dengan kasus distribusi normal bivariat. V.3.3. Fungsi Kumulatif Densitas Empirik (Empirical Cummulative Density Function (ECDF))
R menyediakan fitur dengan fungsi ecdf()untuk mencari fungsi densitas kumulatif empirik (empirical cummulative density function (ecdf)).. Berikut adalah teori yang berkaitan dengan ecdf: • Fungsi distribusi kumulatif suatu populasi dituliskan dalam bentuk
F ( x) = Pr ob( X ≤ x) •
Daerah di bawah fungsi densitas f(x) dari a ke b dituliskan dalam bentuk
F (b) − F (a) = Pr ob(a < X ≤ b) • • • • • •
Fungsi distribusi kumulatif dicari dengan mengestimasi nilai F(x), dimana proporsi nilai data ≤ x Histogram yang ditampilakan adalah histogram kumulatif Dapat menjadi sempurna jika histogramnya hanya memiliki satu observasi per bin (batang) ECDF bersifat unik dan tidak memerlukan binning Sangat baik untuk menunjukkan perbedaan dalam distribusi keseluruhan di antara dua atau tiga grup yang berlapis/ menumpuk. Dapat membaca langsung nilai quantil nya.
V.4. Regresi & ANOVA (Analysis of Variance)
Dalam sub bab ini akan dijelaskan tentang regresi dan ANOVA, yang merupakan salah satu metode dasar statitistik dalam melakukan pengolahan dan analisis data. Disini akan disajikan teori dan praktek regresi dan anova serta penggunaan R untuk analisis tersebut. Pada tahap awal mungkin akan terasa sedikit rumit, namun diharapkan setelah mencoba beberapa contoh akan menjadi lebih mudah. Analisis regresi digunakan untuk menjelaskan atau memodelkan hubungan antara suatu variable tunggal Y, disebut sebagai variable respon, output atau terikat, dan satu atau lebih variable predictor, input, bebas atau penjelasan (explanatory), X1, …, Xp. Apabila p = 1 maka disebut regresi sederhana, sedangkan apabila p > 1 maka disebut regresi berganda atau regresi multivariate. Jika terdapat lebih dari satu variable terikat Y, maka disebut regresi mulrivariate berganda. Variabel respon harus berbentuk kontinu, sedangkan variabel penjelasan dapat berbentuk kontinu, diskrit ataupun kategori. Sebelum kita masuk ke contoh, maka akan digunakan data pima dari paket library ”faraway” yang sudah tersedia di paket
71
CD. Selain dalam CD, library(faraway) dapat didownload di situs www.stat.lsa.umich.edu/˜faraway. Pertama digunakan library faraway dan memanggil data pima tersebut dengan mengetikkan: > library(faraway) > data(pima) > pima pregnant glucose diastolic triceps insulin bmi diabetes age test 1 6 148 72 35 0 33.6 0.627 50 1 2 1 85 66 29 0 26.6 0.351 31 0 3 8 183 64 0 0 23.3 0.672 32 1 4 1 89 66 23 94 28.1 0.167 21 0 5 0 137 40 35 168 43.1 2.288 33 1 ….
Simpulan dari dataset tersebut dapat kita lihat sebagai berikut: > summary(pima) pregnant Min. : 0.000 1st Qu.: 1.000 Median : 3.000 Mean : 3.845 3rd Qu.: 6.000 Max. :17.000 insulin Min. : 0.0 1st Qu.: 0.0 Median : 30.5 Mean : 79.8 3rd Qu.:127.2 Max. :846.0 test Min. :0.0000 1st Qu.:0.0000 Median :0.0000 Mean :0.3490 3rd Qu.:1.0000 Max. :1.0000
glucose Min. : 0.0 1st Qu.: 99.0 Median :117.0 Mean :120.9 3rd Qu.:140.2 Max. :199.0 bmi Min. : 0.00 1st Qu.:27.30 Median :32.00 Mean :31.99 3rd Qu.:36.60 Max. :67.10
diastolic Min. : 0.0 1st Qu.: 62.0 Median : 72.0 Mean : 69.1 3rd Qu.: 80.0 Max. :122.0 diabetes Min. :0.0780 1st Qu.:0.2437 Median :0.3725 Mean :0.4719 3rd Qu.:0.6262 Max. :2.4200
triceps Min. : 0.00 1st Qu.: 0.00 Median :23.00 Mean :20.54 3rd Qu.:32.00 Max. :99.00 age Min. :21.00 1st Qu.:24.00 Median :29.00 Mean :33.24 3rd Qu.:41.00 Max. :81.00
Simpulan data di atas menunjukkan beberapa hal yang tidak wajar yang mungkin disebabkan ketika memasukkan data, misalkan nilai maksimum pregnant = 17 (ini sangat besar, tetapi juga memungkinkan), tekanan darah (diastolic) = 0. Dengan demikian kita perlu sedikit melakukan pembersihan data tersebut, khususnya yang diasumsikan ketika memasukkannya tidak benar. Melihat ketidak wajaran tersebut kita perlu menset nilai dari dataset yang bernilai 0 dengan NA, hal ini untuk menghindari terjadinya pada kesalahan perhitungan. Nilai NA data tidak tersedia untuk nilai yang diasumsikan salah ketika memasukkan, sehingga tidak termasuk dalam perhitungan. Apabila nilai tersebut tetap 0, maka nilai statistik menjadi tidak wajar, karena tidak sesuai dengan kenyataan. Perintah yang digunakan untuk menset nilai tersebut adalah: > > > > >
pima$diastolic[pima$diastolic == 0] <- NA pima$glucose[pima$glucose == 0] <- NA pima$triceps[pima$triceps == 0] <- NA pima$insulin[pima$insulin == 0] <- NA pima$bmi[pima$bmi == 0] <- NA
72
Variabel test adalah bukan kuantitatif melainkan kategori. Variabel tersebut adalah faktor. Tetapi karena merupakan kode numerik, variabel tersebut diperlakukan seperti halnya kuantitatif. Akan lebih baik apabila variabel tersebut adalah faktor sehingga dapat diperlakukan sesuai fungsinya. Hal seperti itu kadang terlupakan dan akibatnya akan terjadi kesalahan proses perhitungan statistik (misalkan proses perhitungan ”rata-rata code pos” (dimana kodenya numerik tetapi bukan data kuantitatif)). Sehingga untuk mendefinisikan variabel test sebagai faktor maka perintah yang digunakan adalah: > pima$test <- factor(pima$test) > summary(pima$test)
0 1 500 268 Sekarang jelas bahwa 500 kasus adalah negative (0) dan 268 positif (1). Sehingga untuk menggunakan deskripsi label yang benar untuk variabel test sebagai kategori maka dimana (0) adalah negatif dan (1) adalah positif, maka dituliskan perintah berikur: > levels(pima$test) <- c("negative","positive") > summary(pima) pregnant Min. : 0.00 1st Qu.: 1.00 Median : 3.00 Mean : 3.85 3rd Qu.: 6.00 Max. :17.00 NA’s : 5 bmi Min. :18.2 1st Qu.:27.5 Median :32.3 Mean :32.5 3rd Qu.:36.6 Max. :67.1 NA’s :11.0
glucose Min. : 44 1st Qu.: 99 Median :117 Mean :122 3rd Qu.:141 Max. :199 NA’s :35 diabetes Min. :0.078 1st Qu.:0.244 Median :0.372 Mean :0.472 3rd Qu.:0.626 Max. :2.420
diastolic Min. : 24.0 1st Qu.: 64.0 Median : 72.0 Mean : 72.4 3rd Qu.: 80.0 Max. :122.0 NA’s :227 age Min. :21.0 1st Qu.:24.0 Median :29.0 Mean :33.2 3rd Qu.:41.0 Max. :81.0
triceps Min. : 7.0 1st Qu.: 22.0 Median : 29.0 Mean : 29.2 3rd Qu.: 36.0 Max. : 99.0 NA’s :374 test negative:500 positive:268
insulin Min. : 14.0 1st Qu.: 76.2 Median :125.0 Mean :155.5 3rd Qu.:190.0 Max. :846.0
Data yang ada menunjukkan beberapa kemungkinan diantaranya regresi test (sebagai respons) terhadap diastolic dan bmi (predictor) melibatkan suatu respon kualitatif. Hal tersebut menunjukkan bahwa analisis regresi memiliki beberapa kemungkinan berdasarkan pada tujuan yang dikehendaki, yaitu: − Prediksi observasi ke depan − Penilaian efek, hubungan antara, variable explanatory terhadap respon. − Deskripsi umum struktur data. Variabel prediktor dalam fungsi adalah berbentuk data kategori/kualitatif. Nama ANOVA sendiri merupakan kependekan dari Analysis of Variance, dimana sering digunakan untuk membagi/mempartisi keseluruhan variansi sebagai respon terhadap setiap faktor-faktor dan error. Pada saat ini variabel prediktor disebut sebagai faktor karena memiliki beberapa tingkatan/level. Sedangkan parameter-paremeter yang terlibat disebut sebagai efek. Kita sekarang akan memperhatikan model dimana parameternya tetap (fixed) tetapi tidak diketahui – disebut model fixed-effect tetapi model randompeffect juga digunakan dimana parameter yang diambil merupakan variabel random. 73
V.4.1. ANOVA Satu Arah V.4.1.1. Model Model untuk ANOVA jenis ini yaitu, diberikan faktor α terjadi pada tingkat i = 0, 1, …, I , dengan sejumlah j observasi per tingkat dimana j = 1, …, Ji. Model yang digunakan: yij = μ + α i + ε ij , i = 1,..., I ; j = 1,..., J i
dengan batasan-batasan parameter sebagai berikut: 1. Set μ = 0 dan variabel dummy I yang berbeda 2. α1 = 0 – berhubungan dengan perlakuan yang kontras 3. ∑ J iα i = 0 , yang mengarah pada penaksir kuadrat
terkecil (least square
i
estimate).
μˆ = y•• , αˆi = yi • − y•• ,
Dimana tanda • (misal pada y•• ) menyatakan indeks dari rerata yang digunakan. Metode terakhir ini paling banyak direkomendasikan untuk penghitungan manualwalaupun akan lebih sulit disajikan dalam bentuk y = Xβ + ε . Faktor kedua pertama di atas lebih mudah diimplementasikan secara komputas. Seperti biasa, beberapa analisis grafis awal diperlukan sebelum melakukan pencocokkan (fitting). Boxplot sisi per sisi merupakan plot yang paling banyak digunakan. Lihat pada persamaan variansi, transformasi, outliers (disini tidak berpengaruh secar relevan karena levarage tidak akan berbeda kecuali disain yang sangat tidak seimbang) V.4.1.2. Penaksiran dan Pengujian
Penaksiran efek menggunakan langsung formula seperti yang sudah dijelaskan sebelumnya atau menggunakan pendekatan least square (karena hasil sama). Pada ANOVA, penaksiran dan pengujian dilakukan pertama kali adalah mencari perbedaan factor levelnya. Kita lakukan perbandingan hipotesis, dimana: • H0 : αi = 0 , berlaku untuk semua i • Ha :paling sedikit satu αi ≠ 0 Contoh kasus, Suatu penelitian ingin mengetahui apakah factor sosioekonomi mempengaruhi IPK mahasiswa. Sosioekonomi dibagi dalam tiga kelompok, yaitu rendah (r), sedang (s) dan tinggi (t). Dalam penelitian tersebut, masing-masing kelompok terdiri dari 7 mahasiswa yang diambil data IPK secara acak. > ipksosek <- edit(data.frame()) #memasukkan data ipksosek >ipksosek 1 2 3 4 5 6
ipk 2.87 2.16 3.14 2.51 1.80 3.01
sosek r r r r r r
74
7 8 9 10 11 12 13 14 15 16 17 18 19 20 21
2.16 3.23 3.45 2.78 3.77 2.97 3.53 3.01 2.25 3.13 2.44 2.54 3.27 2.81 1.36
r s s s s s s s t t t t t t t
Langkah pertama adalah melihat sebaran data dengan melakukan plot, dan digunakan boxplot dan hasilnya seperti Gambar 5.24: > boxplot(ipk˜sosek, data=ipksosek)
Gambar 5.24. Boxplot sebaran IPK mahasiswa berdasarkan sosio ekonomi Disini diharapkan tidak terlihat adanya outlier, kemiringan (skewness) dan variansi yang berbeda. Selanjutnya, dilakukan pencocokan model dengan cara sebagai berikut: Eksperimen tersebut melibatkan suatu factor tunggal yaitu kelas sosioekonomi (sosek) pada tiga level (r,s,t). Misalkan μ r , μ s , μt masing-masing menyatakan ratarata IPK untuk mahasiswa tingkat social ekonomi rendah, sedang dan tinggi. Sehingga akan diujikan hal berikut: H 0 : μ r = μ s = μt melawan alternatif: Ha : paling sedikit dua dari tiga perlakuan berbeda. Maka untuk mendapatkan hal tersebut, maka perlu dibuat suatu model linear yang menggambarkan hubungan ipk dengan sosek sebagai berikut: > g <- lm(ipk˜sosek,ipksosek)
75
> summary(g) Residuals: Min 1Q Median -1.18286 -0.29286 -0.01143
3Q 0.34857
Max 0.72714
Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 2.52143 0.19336 13.040 1.31e-10 *** soseks 0.72714 0.27346 2.659 0.0160 * sosekt 0.02143 0.27346 0.078 0.9384 --Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 0.5116 on 18 degrees of freedom Multiple R-Squared: 0.3372, Adjusted R-squared: 0.2636 F-statistic: 4.579 on 2 and 18 DF, p-value: 0.02468
Karena nilai statistik F = 4.579 melebihi nilai kepercayaan F0.05 = 3.55, maka H0 ditolak dan disimpulkan (pada tingkat kepercayaan α = 0.05) bahwa rata-rata IPK mahasiswa berbeda paling sedikit dua dari tiga tingkat sosio ekonomi. Hal ini juga ditunjukkan oleh α = 0.05 lebih besar dari p-value = 0.02468. Grup r adalah level rujukan dan mempunyai nilai rerata 2.52143, sedangkan grup s dan t berturut-turut memiliki nilai rerata 0.72714 dan 0.02143 lebih besar dari ratarata r. Kemudian berdasarkan factor sosioekonomi tersebut, maka dilakukan pemeriksaan matriks rancangan untuk memahami pengkodean: > model.matrix(g)
Kita dapat mencocokan model tanpa mengintersepsi batasan: > gi <- lm(ipk˜sosek -1, ipksosek) > summary(gi) Residuals: Min 1Q Median -1.18286 -0.29286 -0.01143
3Q 0.34857
Coefficients: Estimate Std. Error t value sosekr 2.5214 0.1934 13.04 soseks 3.2486 0.1934 16.80 sosekt 2.5429 0.1934 13.15 --Signif. codes: 0 ‘***’ 0.001 ‘**’
Max 0.72714
Pr(>|t|) 1.31e-10 *** 1.91e-12 *** 1.14e-10 *** 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.5116 on 18 degrees of freedom Multiple R-Squared: 0.972, Adjusted R-squared: 0.9674 F-statistic: 208.4 on 3 and 18 DF, p-value: 3.662e-14
Kita lihat, dengan nilai p (p-value) = 3.662*10-14 (mendekati 0) maka kita tidak dapat melakukan perbandingan hipotesis sebelumnya.
76
V.4.2. ANOVA Dua Arah
Anggap terdapat dua faktor, α pada level I dan β pada level J. Misalkan jumlah observasi pada level i dari α dan level j dari β dan misalkan observasi tersebut sebagai yi j1, yi j2, …. Suatu layout lengkap A memiliki nij ≥ 1 untuk semua i, j. memiliki model umum yang paling sesuai adalah : yijk = μ + α i + β j + (αβ )ij + ε ijk
Efek interaksi (αβ)ij diinterpretasikan sebagai bagian dari rata-rata respon yang bukan merupak atribut efek tambahan dari αi and βj Sebagai contoh misalkan anda menyukai strawberi dan es krim secara individu, tetapi strawberi yang dicampur dalam eskrim lebih disukai. Hal yang berlawanan misalkan anda menyukai ikan dan es krim tetapi tidak menyukai es krim yang dicampur ikan. (
Berikut ini adalah contoh penggunaan anova dua arah [23]: Perusahaan pengecatan pesawat udara akan melakukan pengecatan pada permukaan alumunium dengan dua cara; penyemprotan (spray) dan dipping. Percobaan perbedaan faktor dilakukan untuk menyelidiki efek dari tipe pengecatan dan metode pengecatan terhadap daya rekat (adhisif) cat. Kekuatan daya rekat diukur, dengan tiga jenis cat (primer) yang berbeda dan dua metode pengecatan (aplikasi). Berikut adalah penulisan dalam R: > paint <- data.frame(adhf = c(4.0,4.5,4.3,5.6,4.9,5.4,3.8,3.7,4.0,5.4,4.9,5.6,5.8,6.1,6.3,5.5 ,5.0,5.0), primer = factor(rep(rep(1:3,rep(3,3)),2)),applic = factor(rep(c("D","S"),c(9,9)))) > paint adhf primer applic 1 4.0 1 D 2 4.5 1 D 3 4.3 1 D 4 5.6 2 D 5 4.9 2 D 6 5.4 2 D 7 3.8 3 D 8 3.7 3 D 9 4.0 3 D 10 5.4 1 S 11 4.9 1 S 12 5.6 1 S 13 5.8 2 S 14 6.1 2 S 15 6.3 2 S 16 5.5 3 S 17 5.0 3 S 18 5.0 3 S
Hasil tersebut adalah untuk anova satu arah, dengan mengabaikan metode pengecatan. Statemen lm() menyesuaikan model linier dan statemen anova() menampilkan hasil dalam bentuk tabel anova. > anova(lm(adhf~primer, data=paint))
77
Analysis of Variance Table Response: adhf
Df Sum Sq Mean Sq F value Pr(>F) 4.5811 2.2906 5.5989 0.01527 * 6.1367 0.4091
primer 2 Residuals 15 --Signif. codes:
0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1
Tampilan hasil di atas adalah anova satu arah untuk metode pengecatan, dengan mengabaikan jenis cat (primer). > anova(lm(adhf~applic, data=paint)) Analysis of Variance Table Response: adhf Df applic 1 Residuals 16 --Signif. codes:
Sum Sq 4.9089 5.8089
Mean Sq 4.9089 0.3631
F value 13.521
Pr(>F) 0.002039 **
0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1
Hasil berikut merupakan anova dua arah untuk metode pengecatan, varibel primer dan interaksi antar keduanya. Sebagai catatan, jumlah kuadrat untuk variable primer dan untuk variable applic bernilai sama dengan hasil perhitungan dalam analisis satu arah. > anova(lm(adhf~primer*applic, data=paint)) Analysis of Variance Table Response: adhf primer applic primer:applic Residuals --Signif. codes:
Df 2 1 2 12
Sum Sq 4.5811 4.9089 0.2411 0.9867
Mean Sq F value 2.2906 27.8581 4.9089 59.7027 0.1206 1.4662 0.0822
Pr(>F) 3.097e-05 *** 5.357e-06 *** 0.2693
0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1
Interaksi variable tidak signifikan (p = 0.2693) oleh karena itu kita dapat melakukan uji untuk efek utama, dimana keduanya memiliki signifikansi yang tinggi. Kesimpulan: Bahwa pemilihan tipe primer dan pemilih metode pengecatan (aplikasi) mempengaruhi kekuatan daya rekat cat, dan selisih antara tiga tipe cat adalah sama untuk kedua cara pengecatan (aplikasi), dan selisih antara dua metode pengecatan adalah sama untuk setiap tipe cat. Untuk uji t contoh pengecatan di atas, dibahas di V.5. sub bab uji t. Perintah berikut digunakan untuk menunjukkan cara perhitungan nilai rerata untuk enam tipe cat yang berbeda: kombinasi variable applic dan cara bagaimana mengatur nilai rerata dalam bentuk matriks untuk memberikan plot interaksi. > split(paint$adhf,paint$applic:paint$primer) $"D:1" [1] 4.0 4.5 4.3 $"D:2" [1] 5.6 4.9 5.4
78
$"D:3" [1] 3.8 3.7 4.0 $"S:1" [1] 5.4 4.9 5.6 $"S:2" [1] 5.8 6.1 6.3 $"S:3" [1] 5.5 5.0 5.0 >sapply (split (paint$adhf, paint$applic:paint$primer), mean) D:1 4.266667
D:2 5.300000
D:3 3.833333
S:1 S:2 S:3 5.300000 6.066667 5.166667
>matrix(sapply(split(paint$adhf, paint$applic:paint$primer),mean),ncol=2) [,1] [,2] [1,] 4.266667 5.300000 [2,] 5.300000 6.066667 [3,] 3.833333 5.166667
Perintah mathplot() akan mem-plot setiap kolom dari matriks pada graph yang sama. > matplot(matrix(sapply(split(paint$adhf, paint$applic:paint$primer),mean),ncol=2), type="l", xlab="Primer",ylab="Daya Rekat")
Gambar 5.25: Plot kolom matrik Dua garis di atas tampak seperti paralel, sesuai dengan kesimpulan sebelumnya bahwa tidak terjadinya suatu interaksi. Oleh karena itu, selisih/perbedaan antara tipe cat (primer) adalah sama untuk setiap metode pengecatan. 79
> matplot(matrix(sapply(split(paint$adhf,paint$primer:paint$applic), mean),ncol=3), type="l", xlab="Aplikasi",ylab="Daya Rekat")
Cara lain, anda dapat juga memplot daya rekat (adhesi) dengan metode pengecatan (aplikasi), dengan garis yang berbeda untuk setiap primer garis tampak parallel (Gambar 5.23).
Gambar 5.26: Matplot daya rekat (adhesi) vs aplikasi V. 5. Analisis (Uji) Statistika t
Berkaitan den gan contoh kasus pengecatan pesawat terbang dengan variable primer (tipe cat) dan variable applic (metode pengecatan; semprot dan dipping) (lihat kasus di bagian anova dua arah), penulisan uji t dan koefisiennya dalam R adalah: > summary(lm(adhf~primer*applic, data=paint)) Call: lm(formula = adhf ~ primer * applic, data = paint) Residuals: Min 1Q -0.40000 -0.16667 Coefficients:
Median 0.03333
3Q 0.21667
Max 0.33333
Estimate Std. Error t value (Intercept) 4.2667 0.1656 25.772 primer2 1.0333 0.2341 4.414 primer3 -0.4333 0.2341 -1.851 applicS 1.0333 0.2341 4.414 primer2:applicS -0.2667 0.3311 -0.805 primer3:applicS 0.3000 0.3311 0.906 --Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05
Pr(>|t|) 7.1e-12 0.000845 0.088949 0.000845 0.436265 0.382736
*** *** . ***
`.' 0.1 ` ' 1
Residual standard error: 0.2867 on 12 degrees of freedom Multiple R-Squared: 0.9079, Adjusted R-squared: 0.8696
80
F-statistic: 23.67 on 5 and 12 DF,
p-value: 7.89e-06
Data di atas menunjukkan nilai uji t dari masing-masing variabel. V.6. R dan Statistika Lanjut V.6.1. Pencocokan distribusi (Fitting Distribution)
Pencocokan distribusi adalah suatu cara untuk menentukan fungsi matematis yang diwakili dalam bentuk variabel-variabel statistika. Masalah yang sering dihadapi oleh statistikawan adalah, jika mempunyai suatu pengamatan/observasi dengan variabel kuantitatif x1, x2, …., xn dan ingin menguji jika pengamatannya itu, menjadi sample dari populasi yang tidak diketahui atau milik dari suatu populasi dengan pdf (fungsi densitas probabilitas) f(x, θ), dimana θ adalah vektor dari parameter yang akan diestimasi dengan data yang ada. Ada empat langkah utama dalam pencocokan distribusi: 1) Pemilihan model: Hipotesis distribusi 2) Estimasi/taksiran parameter 3) Pengukuran Goodnes of Fit 4) Uji Goodness of Fit Empat langkah ini akan dibahas pada sub-bagian berikut ini. V.6.1.1. Pemilihan Model
Langkah pertama dalam pencocokan distribusi adalah memilih model atau fungsi matematis yang mewakili data dengan baik. Semakin baik suatu model maka semakin baik pula keterwakilan data. (Catatan: baik buruknya model, salah satunya dapat diukur dengan error yang dihasilkan). Terkadang untuk pemilihan suatu model atau fungsi, teknik grafik seperti histogram dan lainnya dapat membantu. Namun secara grafik dapat mengakibatkan cenderung subjektif, oleh karena itu ada metode yang didasari pada ekspresi analitik seperti kriteria Pearson’s K. Memecahkan persamaan diferensial tertentu dapat dicapai dengan menggunakan kelompok fungsi yang mewakili distribusi empiris. Kurva dari fungsi tersebut hanya bergantung pada nilai rerata, variabilitas, skewness/kemiringan, dan kurtosis. Dengan menstandarisasi data, tipe kurva yang hanya bergantung pada pengukuran skewness dan kurtosis diformulasikan dalam fungsi berikut: γ 12 (γ 2 + 6) 2 K= 4(4γ 2 − 3γ 12 + 12)(2γ 2 − 3γ 12 ) dimana: n
γ1 =
∑ (x i =1
i
nσ 3
n
γ2 =
∑ (x i =1
− μ )3
i
adalah koefisien skewness Pearson’s
− μ )4
nσ 4
− 3 adalah koefisien kurtosis Pearson’s
81
Berdasarkan nilai K, yang diperoleh dari data yang ada, kita mempunyai suatu fungsi jenis tertentu. Berikut adalah beberapa contoh distribusi kontinu dan distibusi diskrit. Untuk setiap distribusi, terdapat bentuk grafik dan perintah R untuk membentuk grafik tersebut. Untuk data diskrit, dapat melihat distribusi Poisson, dengan fungsi mass probabilitas: f ( x, λ ) = e − λ
λx x!
dimana x = 0, 1, 2, ...
Contoh, penulisan: > x.poi <- rpois(n=200,lambda=2.5) > hist(x.poi,ylab="Frekuensi",xlab="x.poi",main="Distribusi Poisson")
menghasilkan gambar berikut
Gambar 5.27: Histogram distribusi Poisson Untuk data yang kontinu, distribusi yang digunakan antara lain: • Distibusi Normal (Gaussian), dengan fungsi: f ( x, μ , σ ) =
− 1 e 2π σ
1( x − μ ) 2 2σ 2
untuk x ∈ R
Sebagai contoh, dapat dituliskan perintah pembentukan kur va untuk model distribusi Normal (Gaussian): > curve(dnorm(x,m=10,sd=2), from=0,to=20,main="Distribusi Normal")
maka akan menghasilkan tampilan grafik seperti Gambar 5.28 sebagai berikut:
82
Gambar 5.28: Kurva distribusi Normal •
Distribusi Gamma:
f ( x, α , λ ) =
λα α −1 − λx x e untuk x ∈ R + Γ(α )
Sebagai contoh, penulisan: > curve(dgamma(x, scale=1.5, shape=2),from=0, to=15, main="Distribusi Gamma")
Menghasilkan grafik sebagai berikut:
Gambar 5.29: Kurva distribusi Gamma
83
Distribusi Weibull: f ( x,α , β ) = αβ
−α
x
α −1
x [ − ( )α ]
e
β
untuk x ∈ R +
Sebagai contoh, penulisan: > curve(dweibull(x, scale=2.5, shape=1.5),from=0, to=15, main="Distribusi Weibull")
akan menghasilkan grafik seperti Gambar 5.30 sebagai berikut:
Gambar 5.30: Kurva distribusi Weibull Indeks kemiringan dan kurtosis dihitung dengan menggunakan perintah skewness()dan kurtosis() yang terdapat dalam paket fBasics (dapat didownload dari website CRAN). Perintah untuk menghitungnya adalah sebagai berikut: > library(fBasics) ## load paket fBasics > skewness(x.norm) ## kemiringan(skewness) distribusi normal [1] 0.1242952 > kurtosis(x.norm) ## kurtosis distribusi normal [1] 0.01372539 > skewness(x.wei) ## skewness distribusi Weibull [1] 0.7788843 > kurtosis(x.wei) ## kurtosis distribusi Weibull [1] 0.4331281
V.6.1.2. Taksiran parameter
Setelah pemilihan model yang merepresentasikan data, kita perlu melakukan penaksiran parameter dari model yang telah dipilih. Dalam literatur statistik, ada beberapa metode penaksiran parameter namun dalam buku ini memfokuskan pada metode berikut:
84
1) analogi 2) moment 3) maximum likelihood. Metode analogi terdapat dalam penaksiran model parameter yang mengaplikasikan fungsi yang sama untuk data empiris, yaitu: kita estimasi rerata yang tidak diketahui dari populasi normal menggunan rerata sample. > mean.hat<-mean(x.norm) > mean.hat [1] 0.02646597
Metode moment adalalah teknik untuk membentuk penaksir parameter berbasis pada pencocokan moment sample yang berhubungan dengan moment distribusi. Metode ini membandingkan moment terhadap populasi (secara teori). Apabila terdapat metode momen maka keuntungannya adalah terdapat pula bentuk yang sederhana. Kita mendefinisikan moment sample (secara empiris) dengan cara tersebut.
Pendefinisian moment (empiris) suatu sampel adalah sebagai berikut: • Moment sample ke-t sekitar 0 (nol): n
mt = ∑ xit yi
, t = 0,1,2,...
i =1
•
Moment sample ke-t sekitar rerata (mean): n
mt' = ∑ ( xi − μ )t yi
, t = 0,1,2,...
i =1
dimana moment populasi dari tinjauan teori adalah: • Moment populasi ke-t sekitar 0 (nol): α
mt* = ∫ x t f ( x,θ )dx , t = 0,1,2,... β
•
Moment populasi ke-t sekitar rerata: α
m = ∫ ( x − μ )t f ( x,θ )dx , t = 0,1,2,... '* t
β
dimana (β – α) adalah range dimana f(x,θ) didefinisikan, μ adalah rerata dari distribusi, dan yi adalah frekuensi relative empiris, yaitu: kita akan menaksir parameter distribusi gamma menggunakan metode moment dengan mempertimbangkan moment pertama sekitar 0 (maupun rerata) dan moment kedua terhadap rerata (maupun variansi):
α =x λ
α = s2 λ2
85
dimana pada sisi kiri adalah rerata dan variansi distribusi gamma dan sisi kanan adalah rerata sample dan variansi sample terkoreksi. Dengan memecahkan persamaan sebelumnya di atas, didapat: x λˆ = 2 s x2 ˆ α= 2 s Penulisan untuk kasus di atas dengan R adalah sebagai berikut: > x.gam<-rgamma(200,rate=0.5,shape=3.5) # sampling dari distribusi gamma dengan λ = 0.5 dan α = 0.35 (parameter bentuk) > med.gam<-mean(x.gam) # rerata sampel > var.gam<-var(x.gam) # variansi sample > l.est<-med.gam/var.gam # penaksiran lambda (berhubungan dengan rata-rata) > a.est<-((med.gam)^2)/var.gam # penaksiran alfa > l.est [1] 0.579852 > a.est [1] 4.153633
Metode maximum likelihood digunakan dalam inferensi statistic untuk menaksir parameter. Misal kita mempunyai variabel random dengan diketahui pdf f(x,θ) yang menggambarkan karakter kuantitatif dalam populasi. Untuk itu kita harus menaksir vektor konstanta dan parameter tidak diketahui (unknown) θ sesuai data sample: x1, x2, …, xn .
Penaksiran maximum likelihood dimulai dengan model matematika yang diketahui sebagai fungsi likelihood dari data sample. Singkatnya, likelihood dari himpunan data adalah probabilitas untuk memperoleh himpunan data dengan yang sudah diberikan model probabilitas yang dipilih. Ekspresi ini mengandung parameter tidak diketahui (unknown). Nilai-nilai parameter yang memaksimumkan likelihood suatu sampel dikenal sebagai penaksiran maximum likelihood (MLE = Maximum Likelihood Estimate). Definisi fungsi likelihood adalah sebagai berikut: n
L( x1 , x2 ,..., xn ,θ ) = ∏ f ( xi ,θ ) i =1
Selanjutnya kita dapat menggunakan metode analisis matematik (turuan parsial sama dengan nol) bila fungsi likelihoodnya cukup sederhana, namun seringnya kita mengoptimasi L(x1, x2, …, xn, θ) menggunakan metode iterasi. MLE mempunyai beberapa sifat statistik dan beberapa keuntungan. Misalkan, dalam kasus distribusi gamma, fungsi likelihoodnya adalah:
86
n
−λ ∑ x λα α −1 − λx λα n n L( x1 , x2 ,..., xn ,θ ) = ∏ f ( xi ,α , λ ) = ∏ xi e ) (∏ xi )α −1 e =( Γ(α ) i =1 i =1 i =1 Γ (α ) n
n
i
i =1
i
dan fungsi logaritmanya: n
n
i =1
i =1
log( L) = nα log(λ ) − n log(Γ(α )) + (α − 1)∑ log xi − λ ∑ xi
R menyediakan fungsi untuk menghitung MLE dalam dua statement berikut: 1) mle() yang termasuk dalam paket stats4 2) fitdistr() yang termasuk dalam paket MASS mle() mencocokkan (fit) parameter dengan metode maximum likelihood yang menggunakan metode iterasi kalkulus numeric untuk meminimumkan negative loglikelihood ekspresi analitik sebagai argumen dan memberikan nilai penaksiran parameter awal. Penulisan di R untuk kasus distribusi gamma adalah sebagai berikut: > library(stats4) # loading paket stats4 > ll<-function(lambda,alfa) {n<-200 x<-x.gam n*alfa*log(lambda)+n*log(gamma(alfa))-(alfa1)*sum(log(x))+lambda*sum(x)} # fungsi negative log-likelihood > est<-mle(minuslog=ll, start=list(lambda=2,alfa=1)) > summary(est) Maximum likelihood estimation Call: mle(minuslogl = ll, start = list(lambda = 2, alfa = 1)) Coefficients: Estimate Std. Error lambda 0.5290189 0.05430615 alfa 3.6829126 0.35287672 -2 log L: 1044.282
Pertama ditentukan nilai awal parameter penaksir secara acak, tetapi dapat pula digunakan metode moment untuk menentukan nilai awal parameter penaksir tersebut. Perintah mle( ) digunakan untuk menaksir parameter setiap jenis pdf. Untuk hal ini perlu diketahui bahwa ekspresi analitik likelihood perlu diketahui. Dalam paket MASS tersedia fungsi fitdistr() untuk pencocokan maximum likelihood distribusi univariate tanpa informasi tentang likelihood ekspresi analitik. Hal ini cukup dengan menspesifikasikan suatu vector data, tipe pdf (densfun) dan tentu saja daftar (list) nilai awal (start) untuk prosedur iterasi. > library(MASS) # loading package MASS > fitdistr(x.gam,"gamma") # fitting parameter pdf gamma shape 3.68320097 (0.35290545)
rate 0.52910229 (0.05431458)
87
> fitdistr(x.wei,densfun = dweibull, start = list(scale=1,shape=2))# fitting parameter pdf Weibull scale -1.04721828 (0.03814184)
shape 2.04959159 (0.11080258)
> fitdistr(x.norm,"normal") # fitting parameters pdf gaussian mean sd 9.9355373 2.0101691 (0.1421404) (0.1005085)
V.6.1.3. Pengukuran Goodness of Fit
Pengukuran goodness of Fit berguna untuk memasangkan frekuensi secara empiris dengan yang dicocokkan oleh model teoritis. Dalam hal ini terdapat pengukuran absolute dan relative. Diantara pengukuran absolute kita dapat memilih: n
ξ=
∑y i =1
− yi*
i
n n
ξ=
2
∑(y i =1
i
− yi* )
n
dimana yi adalah frekuensi empiric dan yi* adalah nilai yang dicocokkan. Berikut adalah beberapa pengukuran relative: n
δ=
ξ n
∑ yi n
=
∑y
i
i =1
i =1
− yi*
n
∑y
i
i =1
∑ (y n
δ=
ξ
2
2
n
∑ yi n
=
i
i =1
)
2
− yi* / n
n
∑ y /n
i =1
i =1
i
n
2 2
δ=
ξ
2 n
=
∑y i =1
2 i
/n
∑(y i =1
i
n
− yi* ) 2
∑y i =1
2 i
Keterangan: biasanya indeks tersebut merupakan persentase dari rerata yang sesuai. Berikut adalah contoh menggunakan R perhitungan data (berdistribusi Poisson): > > > >
lambda.est<-mean(x.poi) # estimasi lambda parameter tab.os<-table(x.poi) # tabel dengan frekuensi empiris tab.os x.poi
88
0 1 2 3 4 5 6 7 8 21 29 46 53 28 16 4 2 1 > freq.os<-vector() > for(i in 1: length(tab.os)) freq.os[i]<-tab.os[[i]]
#
vektor frekuensi empirik > freq.ex<-(dpois(0:max(x.poi),lambda = lambda.est)*200) # vektor frekuensi yang diharapkan atau fitted (expected) > freq.os [1] 21 29 46 53 28 16 4 2 1 > freq.ex [1] 15.0040080 38.8603808 50.3241931 43.4465534 28.1316433 14.5721912 6.2903292 [8] 2.3274218 0.7535028 > acc<-mean(abs(freq.os-trunc(freq.ex))) # index absolute dari goodness of fit > acc [1] 2.111111 > acc/mean(freq.os)*100 # persentase relative index goodness of fit [1] 17
Suatu teknik grafik untuk menghitung goodness of fit dapat menggambarkan kurva pdf dan histogram dalam satu kesatuan seperti pada gambar 5.31. Adapun penulisannya dalam R adalah sebagai berikut: > > > > > >
h<-hist(x.norm,breaks=15) xhist<-c(min(h$breaks),h$breaks) yhist<-c(0,h$density,0) xfit<-seq(min(x.norm),max(x.norm),length=40) yfit<-dnorm(xfit,mean=mean(x.norm),sd=sd(x.norm)) plot(xhist,yhist, type = "s", ylim = c(0,max(yhist, yfit)), main=”Normal pdf dan histogram”) > lines(xfit,yfit, col=”red”)
89
Gambar 5.31: Kurva dan histogram pdf Normal V.6.1.4. Uji kenormalan
Dalam analisis statitistika, statistikawan sering melakukan pengujian apakah kumpulan data berasal/tidak dari pupulasi Normal, oleh karena itu perlu melakukan uji kenormalan atas data tersebut. Beberapa literatur menyebutkan, ada beberapa uji yang berguna untuk menguji kemiringan (skewness) dan/atau kurtosis dari suatu distribusi berbasiskan pada b3 e b4 (atau gamma3 e gamma4). Salah satu uji yang paling baik dan ampuh untuk menguji kenormalan terutama untuk sample kecil adalah uji Shapiro-Wilk. Dalam uji Shapiro-Wilk, kenormalan diuji dengan mencocokan dua penaksiran variansi alternative, yaitu: penaksir non-parametrik didapat dari kombinasi linear nilai sample terurut dan penaksir parameter biasa. Bobot (ai) dapat dilihat di table statistik: 2
⎛ n ⎞ ⎜ ∑ ai x(i ) ⎟ ⎠ W = ⎝ni =1 ∑ ( xi − x )2 i =1
Fungsi R untuk uji Shapiro-Wilk adalah shapiro.test() yang akan menghitung nilai statitistik W dan p- value. Penulisan perintah di R adalah sebagai berikut: > shapiro.test(x.norm) Shapiro-Wilk normality test data: x.norm W = 0.9938, p-value = 0.5659
90
p-value lebih tinggi dibandingkan tingkat (level) signifikansi yang biasanya digunakan untuk menguji hipotesis statistik. Hipotesis nol akan diterima apabila sample data termasuk dalam distribusi Gaussian. Dalam Ekonometrik, uji Jarque-Bera sering digunakan untuk menguji kenormalan. Uji Jarque-Bera berdasarkan pada pengukuran kemiringan/skewness dan kurtosis dari distribusi dengan mempertimbangakn distribusi asimptotik dari b3 e b4 yang merupakan, dibawah hipotesa nol, suatu chi kuadrat dengan 2 derajat kebebasan. Uji Jarque-Bera di R berada dalam paket tseries (yang harus di download website CRAN) dengan perintah jarque-bera.test() yang akan menghitung nilai statistic, derajat kebebasan, dan p-value. > library(tseries) # loading paket tseries > jarque.bera.test(x.norm) Jarque Bera Test data: x.norm X-squared = 0.539, df = 2, p-value = 0.7638
Suatu teknik pengujian yang diajukan oleh Cucconi (ahli statistik Italia) menguji kenormalan tanpa minimbulkan masalah penaksiran parameter data sample. Misal x1 ≥ x2 ≥, … ≥ xn sample variabel kontinu dan ζ1, ζ2, ζn himpunan bilangan acak normal standard berukuran n, misalkan: n −1
r =ζn
eq=
∑ζ i =1
2 i
n −1
xi − x r + dimana x adalah rerata σˆ n sample dan σˆ adalah akar kuadrat variansi sample. Disini dapat dilihat jika xi berasal dari populasi normal dan yi berdistribusi normal standard maka kita dapat menggunakan uji Kolmogorov-Smirnov untuk menguji hipotesis tersebut. Berikut adalah contoh penulisan di R: kita pertimbangkan transformasi xi: yi = q
> > > > > > >
zz <-rnorm(n=200,m=0,sd=1) # bilangan random sampling dari N(0,1) r <-zz[200] q <-sd(zz[-200]) m <-mean(x.norm) s <-sqrt(var(x.norm)) y <-q*((x.norm-m)/s)+(r/sqrt(200)) ks.test(y,”pnorm”,m=0,sd=1) One-sample Kolmogorov-Smirnov test
data: y D = 0.0298, p-value = 0.9943 alternative hypothesis: two.sided
Paket nortest (harus di download dari situs CRAN) digunakan untuk melakukan uji kenormalan hingga 5 cara yang berbeda, yaitu: 91
1) sf.test()melakukan uji Shapiro-Francia: > library(nortest) #loading paket > sf.test(x.norm) Shapiro-Francia normality test data: x.norm W = 0.9926, p-value = 0.3471
2) ad.test() melakukan uji Anderson-Darling: merupakan modifikasi uji Kolmogorov-Smirnov dan memerlukan distribusi tertentu untuk menghitung nilai kritis. Pada saat ini tabel nilai kritis sudah mencakup distribusi normal, lognormal, exponensial, Weibull, nilai ekstrim tipe I dan logistik. Uji Anderson-Darling berdasarkan pada statistik ini: A2 = − n − S dimana, (2i − 1) [ln F ( xi ) + ln(1 − F ( xn + i +1 ))] , n adalah ukuran sample, dan F(x) adalah n i =1 fungsi distribusi kumulatif (CDF). R hanya menyediakan uji ini untuk cek kenormalan: n
S =∑
> library(nortest) # loading paket > ad.test(x.norm) Anderson-Darling normality test data: x.norm A = 0.4007, p-value = 0.3581
3) cvm.test() melakukan uji Cramer-Von Mises, yang didasari teori statistik: ∞
W2 =
∫ (F ( x) − F ( x) )
2
n
f ( x)dx
−∞
Penulisan uji ini dalam R adalah sebagai berikut: > library(nortest) # loading paket > cvm.test(x.norm) Cramer-von Mises normality test data: x.norm W = 0.0545, p-value = 0.4449
4) lillie.test() melakukan uji Lilliefors: Uji ini juga merupakan modifikasi dari uji Kolmogorov-Smirnov yang tidak bisa digunakan untuk kenormalan apabila rerata dan standar deviasi (yang diambil dari
92
data sample) dari hipotesis distribusi normal tidak diketahui. Uji ini cukup ampuh untuk data sample yang berukuran kecil. Uji Lilliefors mengevaluasi hipotesis bahwa X berdistribusi normal dengan suatu nilai rerata dan variansi berbanding hipotesis alternative bahwa X tidak berdistribusi normal. Uji ini membandingkan distribusi empiric X dengan distribusi normal dengan rerata dan variansi yang sama dengan X. Uji Lilliefors mirip dengan uji Kolmogorov-Smirnov, namun parameter distribusi normalnya ditaksir dari X. Penulisan dalam R adalah sebagai berikut: > library(nortest) # loading paket > lillie.test(x.norm) Lilliefors (Kolmogorov-Smirnov) normality test data: x.norm D = 0.0414, p-value = 0.5509
5) pearson.test() melakukan uji Pearson’s chi-square: Uji ini serupa dengan uji Chi-square χ 2 yang digunakan untuk uji Goodness of fit pada distribusi normal. Berikut adalah penulisan dalam perintah di R untuk uji chisquare: > library(nortest) # loading paket > pearson.test(x.norm) Pearson chi-square normality test data: x.norm P = 10.12, p-value = 0.753
Soal Latihan Bab V
1. Misalkan, tiap akhir tahun perusahaan akan memberikan bonus kepada pegawai berdasarkan penjualan dalam beberapa bulan terkahir dengan nilai (dalam juta) 12 0.4 5 2 50 8 3 1 4 0.25 Pertama buatlah objek penjualan dengan menggunakan fungsi scan() penjualan = scan(),kemudian carilah nilai mean, variansi, standard deviasi, median, minimum, maksimum dan quantile-nya.
93
2. Dengan menggunakan data yang sama, buatlah histogram untuk penjualan tersebut 3. Misalkan, gunakan data set pada R untuk data Orange. Gambarkan boxplot untuk masing-masing variabel (Tree, age, circumference) data tersebut. (petunjuk: > boxplot(nama_variabel, main=”title boxplo”t, horizontal=TRUE) 4. Masukkan data berikut 60 83 71 57 54 47 94 61 25 48 66 84 99 75 89 Buatlah plot stem (stem plot) dan plot daunnya (leaf plot)! 5. Bacalah plot stem dan plot daun berikut, masukkanlah data dan buatlah histogramnya!
6. Pada library (Packages ->Load package ->MASS) MASS, pilih data UScereal yang mengandung informasi tentang sarapan sereal yang populer . Attach data tersebut dengan cara berikut: > > > >
library(‘MASS’) data(‘UScereal’) attach(UScereal) names(UScereal)
Kemudian periksalah hubungan-hubungan berikut dan berik komentar atas hasil tersebut. Anda dapat menggunakan tables, scatterplots, barplots dan lain lain untuk memeriksa. a. Hubungan antara manufacturer dan shelf b. Hubungan antara fat dan vitamins c. Hubungan antara fat dan shelf d. Hubungan antara carbohydrates dan sugars e. Hubungan antara fibre dan manufacturer f. Hubungan antara sodium dan sugars Apakah terdapat hubungan lainnya yang dapat anda prediksi dan periksa? 7. Sebagai upaya meningkatkan kemampuan mengingat mahasiswa, beberapa perguruna tinggi mencoba program pembelajaran sistem block. Misalkan dari 100 mahasiswa dibagi menjadi dua grup masing-masing 50 orang. Satu setengah bagian mencoba program block dan setengah sisanya tidak menggunakan program (nonblock). Diketahui jumlah kehadiran tiap tahun. Akan diuji apakah program block membuat perbedaan pada daya ingat mahasiswa. Datanya sebagai berikut Program Non-Block Block
1 tahun 18 10
2 tahun 15 5
3 tahun 5 7
4 tahun 8 18
> 5 tahun 4 10 94
Lakukan uji hipotesis untuk menentukan apakah ada perbedaan antara dua tipe program tersebut terhadap daya ingat mahasiswa ? 8. Data set InsectSpray pada data built-in di R, berbicara tentang jumlah serangga dalam beberapa area yang diberi perlakukan dengan 6 jenis semprotan (spray) berbeda. Data sudah dalam format untuk dilakukan ANOVA satu arah- yaitu vektor dengan data(count), dan satu factor yang menggambarkan level (spray). Pertama, buatlah boxplot side-by-side untuk melihat apakah ada kesamaan mean. Kemudian lakukan ANOVA satu arah untuk melihat apakah terdapat kesamaan mean ?
95
