ESTIMASI KOEFISIEN KORELASI POLIKORIK MENGGUNAKAN METODE BAYESIAN DENGAN GIBBS SAMPLER Adi Setiawan (
[email protected]) Program Studi Matematika, Fakultas Sains dan Matematika Universitas Kristen Satya Wacana Jl Diponegoro 52-60 Salatiga 50711, Indonesia
Abstract In this paper, it is described how to find the polychoric correlation coefficient by using Bayesian method with Gibbs Sampler. The method is implemented using WinBUGS. Simulation study is done to describe how the method is used. Key-words: twin study, Bayesian method, Gibbs Sampler, polychoric correlation coefficient 1 Pendahuluan Koefisien korelasi polikorik (polychoric correlation coefficient) merupakan suatu ukuran keterkaitan (association) antara dua variable ordinal (Roscino dan Pollice, 2006). Koefisien korelasi tetrakorik merupakan kasus khusus dari korelasi polikorik bila kedua variabel pengamatan ordinal bersifat dikotomi. Pengukuran koefisien ini didasarkan pada anggapan bahwa dua variabel laten yang berdistribusi normal bivariat menghasilkan pasangan skor ordinal. Pasangan skor ordinal ini dapat ditemukan pada sifat biologis (trait) pada pasangan kembar dalam studi pasangan kembar (twin study). Dalam makalah ini, disajikan bagaimana metode Bayesian dengan Gibbs Sampler digunakan untuk mengestimasi koefisien korelasi polikorik berdasarkan data simulasi yang nantinya dapat digunakan dalam studi pasangan kembar. 2 Dasar Teori Misalkan Y1 dan Y2 adalah ukuran dari suatu trait polikotomi pada 2 anggota pasangan kembar Kita menganggap bahwa vektor (Y1,Y2)t tergantung pada variabel laten (X1, X2)t dan suatu batas (threshold) b1 dan b2 melalui persamaan Yi = 1 jika Xi b1, = 2 jika b1 < Xi b2, = 3 jika Xi > b1 untuk i=1,2 Diasumsikan bahwa (X1, X2)t mempunyai distribusi normal bivariat dan dapat didekomposisi menjadi : X 1 A1 E1 (2) X 2 A2 E2 dengan 0 2 2 A1 , ~ N 2 , 2 0 2 A2
0 2 0 E1 . ~ N 2 , 0 0 2 E2 t t Dalam hal ini (A1,A2) dan (E1,E2) saling bebas dan [-1,1] Matriks (X1,X2)t didekomposisi menjadi dua bagian yaitu (A1,A2)t dan (E1,E2)t dengan adalah bagian dari (X1,X2)t sehingga Corr(A1,A2) = dan (E1,E2)t adalah bagian dari (X1,X2)t sehingga E1 dan E2 saling bebas Koefisien korelasi polikorik didefinisikan sebagai koefisien korelasi antara X1 dan X2 yaitu Cov( X 1 , X 2 ) 2 = Corr(X1,X2) = 2 2 V ( X1) V ( X 2 ) Probabilitas bersyarat bahwa Y1 = 1 diberikan A1 dan A2 adalah P( Y1 1| A1 , A2 ) P( X 1 b1 | A1 , A2 ) P( A1 E1 b1 | A1 , A2 ) P( E1 b1 A1 | A1 , A2 ) b a 1 1 2
,
probabilitas bersyarat bahwa Y1 = 2 diberikan A1 dan A2 adalah P( Y1 2 | A1 , A2 ) P( b1 X 1 b1 | A1 , A2 ) P(b1 A1 E1 b2 | A1 , A2 ) P( b1 A1 E1 b2 A1 | A1 , A2 ) b a1 b a 2 1 1 2 2
,
dan probabilitas bersyarat bahwa Y1 = 3 diberikan A1 dan A2 adalah P( Y1 3| A1 , A2 ) P( X 1 b2 | A1 , A2 ) P( A1 E1 b2 | A1 , A2 ) P( E1 b2 A1 | A1 , A2 ) b a 1 2 1 2
Jika diberikan A1 dan A2 maka variabel Y1 dan Y2 saling bebas Probabilitas bersyarat dari Y2 jika diberikan A1 dan A2 dapat ditentukan dengan cara yang sama Batas bj dapat distandardisasi menjadi bj bj bj (3) V ( X1 ) 2 2 untuk j = 1,2 Misalkan pemetaan (y1, y2, a1, a2) p(y1, y2, a1, a2) adalah fungsi densitas dari vektor (Y1, Y2, A1, A2 ) dan yj p(yj | a1, a2) adalah densitas bersyarat dari Yj diberikan A1 dan A2 untuk j =1,2 Fungsi likelihood untuk vektor (Y1, Y2, A1, A2) dalam n pasangan kembar adalah n
L p ( yi1 , yi 2 , ai1 , ai 2 ) i 1 n
p ( yi1 , yi 2 | ai1 , ai 2 ) f (ai1 , ai 2 ) i 1
n
q ( yi1 , ai1 ) q ( yi 2 , ai 2 ) g (ai1 | ai 2 ) h(ai 2 ) i 1
dengan f adalah densitas gabungan dari (A1,A2), g adalah densitas bersyarat dari A1 diberikan A2 dan h adalah densitas marginal dari A2 yang masing-masing diberikan oleh 1
f ( a i1 , a i 2 )
2
g ( a i1 | ai 2 )
h(ai 2 )
2
1
2
1 2 2 exp 2 a i1 2 a i1 a i 2 a i 2 2 2 1
( a ai 2 ) 2 , exp i1 2 2 2 (1 2 ) 2 2 1 1
a 2 exp i 2 2 , 2 2 2 1
dan 1 y
b a q( y, a ) 2
b a 1 2
y
(b a) 1 1 y 1 (b a) 1 y dengan 1= 12 Dalam hal ini, (A1,A2) tidak teramati (unobservable) Fungsi likelihood yang diberikan dapat dianggap sebagai likelihood lengkap (full likelihood) dengan likelihood yang sebenarnya dimarginalisasi atas missing data Untuk mengkonstruksikan Gibbs sampler, diperlukan fungsi likelihood lengkap Fungsi likelihood ini akan sebanding dengan n
L [ 2 ] n / 2 [1 2 ] n / 2 exp x q ( yi1 , ai1 ) q( yi 2 , ai 2 ) i 1
untuk x
n i 1
(ai1 ai 2 ) 2
n 2 i 1 i 2 2
a
. 2(1 ) 2 Dengan kata lain, fungsi likelihood untuk n pasangan kembar akan sebanding dengan L 2
n / 2
2
2
n
[1 2 ]n / 2 exp z q ( yi1 , ai1 ) q ( yi 2 , ai 2 ) i 1
untuk
2 i1 (ai1 ai 2 ) 2 2 i1 ai 2 z 2(1 2 ) 2 n
n
2
dan 2= 12 Prior konjugat untuk parameter 2= 12 dan 1= 12 dipilih dari keluarga distribusi gamma sehingga fungsi densitas priornya adalah p1 p p11 1 (1 ) 2 1 exp( p21 ) , ( p1 ) p3
p p 3 1 2 ( 2 ) 4 2 exp( p4 2 ) , ( p3 ) distribusi prior untuk parameter adalah 1 3 ( ) I[ 1,1] ( ) , 2
dan distribusi konjugat prior untuk parameter b1 dan b2 dari keluarga distribusi normal sehingga fungsi densitas prior : 1 (b1 p5 ) 2 4 (b1 ) exp( ) 2 p6 2 p6
5 (b2 )
1 (b p7 ) 2 exp( 2 ) 2 p8 2 p8
Dalam hal ini p1, p2, , p8 adalah parameter dipilih yang sesuai Berdasarkan pada anggapan bahwa parameter saling bebas maka fungsi densitas bersamanya adalah (1 , 2 , , b1 , b2 ) 1 (1 ) 2 ( 2 ) 3 ( ) 4 (b1 ) 5 (b2 ) Akibatnya fungsi densitasnya sebanding dengan (1 , 2 , , b1 , b2 ) (1 , 2 , , b1 , b2 ) L . Fungsi densitas posterior bersama (1 , 2 , , b1 , b2 | data ) memenuhi ( p1 n 1)
(1 , 2 , , b1 , b2 | data ) 1
(1 2 )
n 2
exp[ w1 ] w2 ,
dengan w1 p2 1 u 2
(b1 p5 ) 2 (b2 p7 ) 2 , 2 p6 2 p8
n
w2 q ( yi1 , ai1 , ci ) q ( yi 2 , ai 2 , c) i 1
dan u p4
n 1 1 n 2 2 ( a a ) ai 2 i 1 i 2 2 2(1 ) i 1 2 i 1
p4 v untuk
1 n n n 2 2 a 2 i 1 ai1ai 2 i 1 ai 2 2 i 1 i1 2(1 ) Berdasarkan pada fungsi densitas bersama (1 , 2 , , b1 , b2 | data ) , distribusi bersyarat penuh (full conditional distribution) untuk masing-masing parameter dapat dinyatakan sebagai berikut : v
( p1 n 1)
(1 | yang lain) 1
n
exp[ p21 ] q ( yi1 , ai1 ) q ( yi 2 , ai 2 ) , i 1
( p 3 n 1)
( 2 | yang lain) 2 exp[u 2 ] , 2 n / 2 ( | yang lain) (1 ) exp[v 2 ] I [ 1,1] ( ) , (b1 p5 ) 2 (b1 | yang lain) exp w2 , 2 p6 (b p7 ) 2 (b2 | yang lain) exp 2 w2 2 p8 Untuk variabel laten, distribusi bersyarat penuh dapat ditentukan dengan (a ai 2 ) 2 2 (ai1 | yang lain) exp i1 q( yi1 , ai1 ) 2(1 2 ) dan
(a ai1 ) 2 2 (ai 2 | yang lain) exp i 2 q( yi 2 , ai 2 ) 2(1 2 ) Untuk mengkonstruksikan Gibbs Sampler digunakan algoritma sebagai berikut : 1 Inisialisasi parameter [2]0, 0 , [b1]0, [b2]0, [1]0, [a11]0, [a12]0, , [an1]0, [an2]0 dan diset j=1 2 Dibangkitkan [1]j 1 (1 | yang lain) dengan yang lain berarti parameter yang lain 3 Dibangkitkan j (| yang lain) 3 Dibangkitkan [2]j 2 (2 | yang lain) 4 Dibangkitkan [b1]j b1 (b1 | yang lain) 6 Dibangkitkan [b2]j b2 (b2 | yang lain) 7 Dibangkitkan [ai1]j ai1 (ai1 | yang lain) 8 Dibangkitkan [ai2]j ai2 (ai2 | yang lain) 9 Langkah 2 sampai 8 untuk j = 2, 3, sampai rantai Markov (Markov chain) hasil dari Gibbs Sampler menjadi konvergen Distribusi bersyarat penuh dari parameter 2 merupakan distribusi ( p3+n, u) yang merupakan distribusi standard sehingga sampling dari parameter 2 mudah diimplementasikan Distribusi bersyarat penuh yang lain tidak ada yang merupakan anggota keluarga standard Untuk menyampelnya digunakan algoritma Metropolis-Hasting dengan densitas yang diusulkan adalah densitas eksponensial dengan mean 1 yaitu y 1 p( y |1 ) exp , y 0 1 1 Distribusi bersyarat penuh untuk juga bukan merupakan distribusi non-standard Untuk memperoleh sampelnya kita menggunakan algoritma Metropolis-Hasting dengan densitas yang diusulkan 1 q( y | ) 2 dengan min{ - , -1} < y < max{ + , 1} dan sebagai contoh diambil = 03 Lebih jauh, distribusi bersyarat penuh dari b1 juga merupakan distribusi non standard Distribusinya dapat didekati dengan menggunakan algoritma Metropolis-Hasting dengan distribusi proposal N(b1,1) Parameter b2, ai1 dan ai2 untuk i = 1, 2, , n dapat dikerjakan dengan cara yang sama 3 Studi Simulasi, Hasil dan Pembahasan Paket program R digunakan untuk membangkitkan data kategorikal pada sampel pasangan kembar yang dianggap tergantung pada variabel laten yang berdistribusi normal bivariat baku dengan koefisien korelasi . Dalam simulasi ini, dipilih untuk menggunakan 600 pasangan kembar dan menggunakan input koefisien korelasi polikorik = 0,2; 0,5; 0,8 sebagai input yang digunakan untuk menjelaskan penggunaan metode tersebut. Tabel 1 merupakan contoh dari hasil simulasi tersebut bila digunakan = 0,2. Dari tabel tersebut, terlihat bahwa dari 600 pasang kembar,
terdapat 27 pasang dengan kembar 1 mempunyai kategori 1 sedangkan kembar 2 mempunyai kategori 1, 59 pasang dengan kembar 1 mempunyai kategori 1 sedangkan kembar 2 mempunyai kategori 2 dan seterusnya. Tabel 1. Tabel kontingensi dari banyaknya kembar 1 dan pasangannya kembar 2 yang berstatus kategori 1, kategori 2 dan kategori 3
Kembar 1 Kategori 1 Kategori 2 Kategori 3
Kategori 1 27 57 9
Kembar 2 Kategori 2 59 269 86
Kategori 3 8 66 8
Bila simulasi dilakukan n = 5 kali maka akan diperoleh hasil lengkap seperti dinyatakan pada Tabel 2. Untuk mengestimasi besarnya koefisien polikorik dengan metode yang telah dijelaskan di atas diimplementasikan dalam WinBUGS versi 14 Prior yang digunakan untuk parameter 12 dan 12 adalah prior (1,1) Berdasarkan pada prior ini, variansi dari variabel laten X1 yang dinyatakan dengan V(X1) = 2 + 2 akan memberikan probabilitas yang tinggi pada interval (0,10) Pada prior dari b1 digunakan distribusi N(0,1) sehingga dalam pandangan persamaan (3), parameter b1 memberikan probabilitas yang tinggi pada interval ( -310, 310) sehingga prior distribusi N( 0, 10) akan merupakan pemilihan yang beralasan untuk parameter b1 Dalam hal ini digunakan median dari rantai Markov dalam MCMC untuk mengestimasi b1, b2 dan Nilai dalam tanda kurung memberikan estimasi interval kredibel (credible interval) 95 % untuk hasil estimasi Bayesian Gambar 1 dan Gambar 2 masing-masing memberikan plot MCMC dan estimasi density (kernel density estimation) untuk parameter-parameter yang diperlukan Berdasarkan pada Tabel 2, terlihat bahwa metode Bayesian memberikan estimasi yang relatif memuaskan karena sesuai dengan parameter yang digunakan untuk membangkitkan data kategorikal pasangan kembar. Tabel 2. Data hasil simulasi dari trait kategorikal pada pasangan kembar. No 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
0,2 0,2 0,2 0,2 0,2 0,5 0,5 0,5 0,5 0,5 0,8 0,8 0,8 0,8 0,8
(1,1) 27 23 20 21 25 49 41 40 40 42 47 41 59 67 56
(1,2) 59 62 59 63 56 56 63 46 39 49 38 29 37 38 36
(1,3) 8 9 5 8 9 3 3 1 1 3 0 0 0 0 0
(2,1) 57 58 56 58 52 59 58 53 59 62 33 39 47 42 36
(2,2) 269 292 302 283 296 296 290 310 326 304 352 353 339 338 332
(2,3) 66 64 62 66 72 50 41 60 42 54 28 36 30 39 52
(3,1) 9 7 9 7 7 1 3 3 2 2 0 0 0 0 0
(3,2) 86 58 64 67 66 52 56 41 41 54 27 39 34 29 35
(3,3) 8 9 5 8 9 3 3 1 1 3 0 0 0 0 0
0 2 4 6 8
-1.00 -1.20
0
2000
4000
6000
8000
10000
-1.2
-1.1
-1.0
-0.9
b1
0.8
1.0
0 2 4 6 8
b1
0
2000
4000
6000
8000
10000
0.8
0.9
1.1
b2
0
2
4
6
0.05 0.20 0.35
b2
1.0
0
2000
4000
6000
8000
10000
Tau
0.1
0.2
0.3
0.4
Tau
Gambar 1 Plot MCMC dan estimasi densitas dengan ukuran sampel 10000 untuk batas b1, b2 dan Tabel 3. Hasil estimasi korelasi polikorik dengan menggunakan metode Bayesian. No 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
b1 -1,01 (-1,10 , -0,92) -0,95 (-1,05 , -0,87) -1,03 (-1,12 , -0,94) -1,04 (-1,14, -0,95) -1,06 (-1,15 , -0,97) -0,91 (-1,01 , -0,81) -0,93 (-1,03 , -0,84) -1,03 (-1,12 , -0,93) -1,02 (-1,13 , -0,92) -0,96 (-1,07 , -0,87) -1,09 ( -1,19, -0,98) -1,15 (-1,26, -1,02) -0,96 (-1,06 , -0,86) -0,93 (-1,03, -0,84 ) -1,02 (-1,13 , -0,92)
b2 0,94 (0,85 , 1,03) 0,98 (0,90 , 1,07) 0,99 ( 0,90 , 1,09) 0,96 (0,87 , 1,05) 1,01 ( 0,91 , 1,10) 1,06 (0,96 , 1,16) 0,99 (0,89 , 1,09) 0,98 (0,88 , 1,07) 1,02 (0,91 , 1,12) 1,06 (0,96 , 1,16) 0,96 (0,85 , 1,06) 0,97 (0,87 , 1,08) 1,07 (0,96 , 1,17) 1,10 (1,01 , 1,19) 0,99 (0,89 , 1,10)
0,18 (0,08 , 0,29) 0,23 (0,14 , 0,33) 0,23 (0,13 , 0,33) 0,22 (0,13 , 0,32) 0,17 (0,06 , 0,30) 0,53 (0,44 , 0,60) 0,52 (0,43 , 0,60) 0,56 (0,49, 0,63) 0,63 (0,55, 0,69) 0,48 (0,39 , 0,56) 0,84 (0,78 , 0,88) 0,78 (0,72 , 0,84) 0,79 (0,72 , 0,84) 0,80 (1,01 , 1,19) 0,76 (0,69 , 0,82)
4 Kesimpulan dan Saran Dalam makalah ini telah dijelaskan bagaimana menggunakan metode Bayesian dengan Gibbs Sampler untuk menentukan koefisien korelasi polikorik. Dengan menggunakan data hasil simulasi yang dibangkitkan menggunakan input
yang dipilih maka nilai-nilai input tersebut dapat diestimasi ulang dengan menggunakan metode Bayesian. Beberapa penelitian yang relatif baru terkait dengan penggunaan metode Bayesian dalam studi twin dapat dilihat pada makalah Eaves dan Erkanli (2003), van den Berg et al. (2006) dan Setiawan (2008). Penelitian lebih lanjut dapat dilakukan untuk data-data real pada studi twin atau family study. 5 Daftar Pustaka [1] Berg, S. M. van den, Setiawan, A., Bartels, M., Polderman, T.J.C., van der Vaart, A.W., Boomsma, D.I., (2006), Individual Differences in Puberty Onset in Girls : Bayesian Estimation of Heritabilities and Genetic Correlations, Behavior Genetics, 36 (2) : 261-270. [2] Cowles, M. K., (2004), Review of WinBUGS 1.4, Am. Stat. 58:330-336. [3] Eaves, L. J., dan Erkanli, A., (2003) Markov Chain Monte Carlo approaches to analysis of genetic and environmental components of human developmental change and G × E interaction. Behav. Genet. 33:279-299. [4] Roscino, A. dan A. Pollice, (2006) A Generalization of the Polychoric Correlation Coeffiecient, Data Analysis, Classification and the Forward Search, Springer, Berlin. [5] Setiawan, A., (2008) Estimasi Bayesian Untuk Penentuan Besarnya Pengaruh Genetik Terhadap Sifat Fenotip dan Studi Simulasinya, Prosiding Seminar Nasional Matematika dan Pendidikan Matematika ISBN : 978-979-16353-1-8.
