MODEL REGRESI LATEN PADA EFEK PLASEBO
DIANA PURWANDARI
DEPARTEMEN MATEMATIKA FAKULTAS MATEMATIKA DAN ILMU PENGETAHUAN ALAM INSTITUT PERTANIAN BOGOR BOGOR 2011
ABSTRAK DIANA PURWANDARI. Model Regresi Laten pada Efek Plasebo. Dibimbing oleh ENDAR HASAFAH NUGRAHANI dan NGAKAN KOMANG KUTHA ARDANA. Analisis regresi mempelajari bentuk hubungan antara peubah respons dengan satu atau lebih peubah prediktor. Analisis akan bertambah kompleks ketika melibatkan prediktor laten, yaitu peubah yang tidak teramati. Model regresi dengan prediktor laten disebut model regresi laten. Model regresi laten memainkan peran penting dalam pemodelan data dengan peubah laten, misalnya efek plasebo pada penyembuhan depresi. Pada model regresi laten ini, prediktor laten diasumsikan kontinu dan diduga menggunakan distribusi beta. Parameter model diduga dengan algoritma EM (Expectation-Maximization). Implementasi model dilakukan dengan menggunakan software R 2.13.1. Hasil analisis data suatu penelitian efek plasebo pada penyembuhan depresi menunjukkan tingkat kesesuaian yang tinggi. Kata kunci: regresi laten, algoritma EM, distribusi beta, efek plasebo.
ABSTRACT DIANA PURWANDARI. Latent Regression Model on Placebo Effect. Under supervision of ENDAR HASAFAH NUGRAHANI and NGAKAN KOMANG KUTHA ARDANA. Regression analysis models the relationship between response variables with one or more predictor variables. The complexity of the analysis will increase, when it involves latent predictors, which are unobserved. Regression model with latent predictor is called latent regression model. Latent regression model has an important role in data modelling with latent variables, such as placebo effect in healing of depression. In this latent regression model, the latent predictor is assumed to be continue and is estimated by using beta distribution. The parameters of the model are estimated by EM (Expectation-Maximization) algorithm. This algorithm is implemented using R 2.13.1 software. The result of data analysis of a placebo effect research on the healing of depression shows a high level of concordance. Keywords: latent regression, EM algorithm, beta distribution, placebo effect.
MODEL REGRESI LATEN PADA EFEK PLASEBO
DIANA PURWANDARI
Skripsi sebagai salah satu syarat untuk memperoleh gelar Sarjana Sains pada Departemen Matematika
DEPARTEMEN MATEMATIKA FAKULTAS MATEMATIKA DAN ILMU PENGETAHUAN ALAM INSTITUT PERTANIAN BOGOR BOGOR 2011
Judul : Model Regresi Laten Pada Efek Plasebo Nama : Diana Purwandari NIM : G54070011
Disetujui Pembimbing I
Pembimbing II
Dr. Ir. Endar H. Nugrahani, MS. NIP. 19631228 198903 2 001
Ir. Ngakan Komang Kutha Ardana, M.Sc. NIP. 19640823 198903 1 001
Diketahui Ketua Departemen Matematika
Dr. Dra. Berlian Setiawaty, MS. NIP. 19650505 198903 2 004
Tanggal Lulus : ………………………………
KATA PENGANTAR Puji dan syukur penulis panjatkan kepada Allah SWT atas segala rahmat dan karunia-Nya serta shalawat dan salam kepada Nabi Muhammad SAW sehingga karya ilmiah ini berhasil diselesaikan. Penyusunan karya ilmiah ini juga tidak lepas dari peranan berbagai pihak. Untuk itu penulis mengucapkan terima kasih yang sebesar-besarnya kepada: 1. Orang tua tercinta: Ateng Rachmat dan Ida Nurdiana yang telah mencurahkan kasih sayangnya, doa, dukungan, kesabaran, kepercayaan, adik Riandi Angga Permana dan Ghina Shoda Nabilah, Om Risdiana serta keluarga besar baik dari papah maupun dari mamah. 2. Dr. Ir. Endar Hasafah Nugrahani, M.S selaku dosen pembimbing I yang telah membimbing, memberikan ilmu, kesabaran, motivasi, dan bantuannya selama penulisan skripsi ini. 3. Ir. Ngakan Komang Kutha Ardana, M.Sc selaku dosen pembimbing II yang telah membimbing, memberikan ilmu, saran, motivasi, dan bantuannya selama penulisan skripsi ini. 4. Dr. Ir. Budi Suharjo, MS selaku dosen penguji yang telah memberikan saran yang terkait dengan skripsi. 5. Staf Departemen Matematika: Pak Yono, Bu Susi, Mas Heri, Pak Bono, Bu Ade, Mas Deni, dan lainnya yang telah memberikan banyak bantuan dan dukungan. 6. Teman-Teman Matematika 44 yang mendukung proses penyusunan karya ilmiah ini: Lukmanul Hakim, Dwi Tanty, Herlan, Wahyu Sudrajat, Yanti AA. 7. Kakak kelas angkatan 41, 42, dan 43 yang selalu menjadi panutan baik. 8. Teman-teman Matematika angkatan 44: Wahyu, Ayum, Iam, Fikri, Wenti, Ririh, Sri, Fajar, Mutia, Rachma, Ayung, Iful, Cita, Tanty, Arina, Deva, Yuyun, Lingga, Masayu, Ruhiyat, Yogie, Lugina, Yanti, Selvie, Nurul, Pepi, Devi, Istiti, Iresa, Sari, Anis, Aqil, Lilis, Imam, Aswin, Eka, Aze, Ali, Vianey, Nadiroh, Nurus, Na’im, Dhika, Ima, Dora, Atik, Nunuy, Yuli, Fani, Phunny, Dian, Rofi, Della, Tyas, Denda, Pandi, Rizqy, Indin, Sholih, Siska, Lili, Tita, Lina, Endro, Lukman, Puying, Tendhy, Ikhsan, Chopa, dan Zae. 9. Adik-adik Matematika angkatan 45 dan 46 yang selalu memberikan dukungan. 10. Sahabatku: Yanti Anjarwati Abbas, Masayu Nur Dzikriyana, Indin Fabrina F, Maryam, Chichi Ryzki, Noe dan Norita yang memberikan semangat, motivasi dan saran-saran. 11. Teman-teman HIMALAYA: Dita, Dwi, Dede Hermanudin, kakak Aris, kakak Arif dan lainnya 12. Teman-teman kosan ceriwis: Ari dan Fahri atas motivasi dan kebersamaannya. 13. Semua pihak yang telah membantu sehingga bisa terselesaikan karya ilmiah ini. Semoga karya ilmiah ini dapat bermanfaat bagi dunia ilmu pengetahuan khususnya matematika dan menjadi inspirasi bagi penelitian-penelitian selanjutnya.
Bogor, November 2011
Diana Purwandari
RIWAYAT HIDUP Penulis dilahirkan di Tasikmalaya pada tanggal 10 Maret 1990 dari pasangan Ateng Rachmat dan Ida Nurdiana. Penulis merupakan anak pertama dari tiga bersaudara. Pendidikan yang telah ditempuh oleh penulis antara lain SDN Kawalu 2 tahun 1995-2001, SMPN 2 Tasikmalaya tahun 2001-2004, SMAN 1 Tasikmalaya tahun 2004-2007. Penulis diterima di Institut Pertanian Bogor melalui jalur USMI (Undangan Seleksi Mahasiswa IPB) pada tahun 2007. Penulis memilih mayor Matematika dengan minor Statistika Terapan, Departemen Matematika, Fakultas Matematika dan Ilmu Pengetahuan Alam. Selama mengikuti perkuliahan, penulis menjadi asisten mata kuliah Kalkulus II (S1) pada semester ganjil tahun akademik 2009-2010, dan asisten praktikum mata kuliah fisika dasar I (S1) pada semester ganjil dan genap tahun akademik 2009-2010. Pada tahun 2007, penulis memperoleh juara 1 lomba baca puisi tingkat asrama TPB IPB dan juara 2 lomba baca puisi islami IPB. Pada tahun 2010, penulis memperoleh juara 2 lomba baca puisi IPB Art Contest dan juara 3 musikalisasi drama FMIPA IPB. Pada tahun 2011, penulis memperoleh juara 1 lomba baca puisi IPB Art Contest dan juara 1 lomba baca puisi IPB NEO EKSMUS. Penulis aktif di berbagai kegiatan kemahasiwaan. Penulis pernah memegang amanah sebagai Wakil Bendahara 1 di Gugusan Mahasiswa Matematika (GUMATIKA) tahun 2008-2009. Penulis juga pernah aktif di Himpunan Mahasiswa Tasikmalaya (HIMALAYA) dan Teater Rumput Fakultas Kehutanan.
DAFTAR ISI Halaman DAFTAR GAMBAR ............................................................................................................ viii DAFTAR LAMPIRAN ......................................................................................................... viii I
PENDAHULUAN ........................................................................................................ 1.1 Latar Belakang ..................................................................................................... 1.2 Tujuan .................................................................................................................. 1.3 Ruang Lingkup .....................................................................................................
1 1 1 1
II
TINJAUAN PUSTAKA ............................................................................................... 2.1 Analisis Regresi .................................................................................................... 2.2 Peubah Acak dan Distribusi .................................................................................. 2.3 Distribusi Bernoulli .............................................................................................. 2.4 Family Beta dan Distribusi Beta ............................................................................ 2.5 Matriks ................................................................................................................. 2.6 Metode Maximum Likelihood ................................................................................ 2.7 Algoritma K-Means ..............................................................................................
2 2 3 3 3 3 4 4
III
METODE PENELITIAN............................................................................................... 3.1 Studi Literatur ...................................................................................................... 3.2 Sumber Data ......................................................................................................... 3.3 Analisis dan Pemrograman ....................................................................................
4 4 5 5
IV
PEMBAHASAN ........................................................................................................... 4.1 Efek Plasebo ......................................................................................................... 4.2 Model Regresi Laten ............................................................................................. 4.3 Metode Pendugaan Parameter ............................................................................... 4.4 Aplikasi pada Efek Plasebo ...................................................................................
5 5 5 6 7
V
SIMPULAN DAN SARAN .......................................................................................... 5.1 Simpulan .............................................................................................................. 5.2 Saran ....................................................................................................................
9 9 9
DAFTAR PUSTAKA ...........................................................................................................
9
LAMPIRAN .........................................................................................................................
10
vii
DAFTAR GAMBAR 1
Halaman Plot kepadatan beta dengan a = 0.5 dan b = 0.3 untuk variabel laten x ............................ 7
2
Histogram data pasien depresi rawat jalan .....................................................................
7
3
Histogram penurunan depresi ........................................................................................
8
4
Grafik estimasi model regresi laten ................................................................................
8
5
Plot kuantil-kuantil model regresi laten untuk data depresi .............................................
8
DAFTAR LAMPIRAN Halaman 1
Pembuktian Teorema Estimasi 𝜷 ...................................................................................
11
2
Penjelasan Model Regresi Linear Multivariat ................................................................
13
3
Pembuktian Persamaan (24) ...........................................................................................
14
4
Pembuktian Persamaan (26) ..........................................................................................
15
5
Pembuktian Teorema Trace ...........................................................................................
16
6
Pembuktian Pendugaan Parameter 𝜷 dengan Algoritma EM (Expectation-Maximization)
18
7
Kode R untuk Gambar 1 ................................................................................................
22
8
Kode R sebagai Pendefinisian Awal (Tarpey & Petkova) ................................................
23
9
Kode R untuk Algoritma K-Means, Maximum Likelihood dan Algoritma EM (Tarpey & Petkova) ........................................................................................................................
25
10 Kode R untuk Plot Kuantil-Kuantil Data Depresi ...........................................................
28
viii
I PENDAHULUAN 1.1 Latar Belakang Penelitian Analisis regresi mempelajari bentuk hubungan antara peubah respons (y) dengan satu atau lebih peubah prediktor (x) terutama untuk menelusuri pola hubungan yang modelnya belum diketahui. Dalam analisis regresi dipelajari bagaimana peubah tersebut berhubungan dan dinyatakan dalam persamaan matematika. Untuk melihat hubungan antara peubah respons dan peubah prediktor secara simultan dapat menggunakan analisis regresi linear. Analisis regresi linear dibedakan menjadi dua, yaitu analisis regresi linear sederhana dan analisis regresi linear berganda. Perbedaan kedua analisis tersebut terletak pada banyaknya peubah prediktor. Analisis regresi linear sederhana terdapat satu peubah prediktor, sedangkan analisis regresi linear berganda terdapat lebih dari satu peubah prediktor (Draper & Smith 1992). Analisis akan bertambah kompleks ketika melibatkan peubah laten, yaitu peubah yang tidak teramati (unobserved variable). Misalkan nilai kuantitas pada peubah laten diperoleh melalui prosedur estimasi (perkiraan). Sebuah model regresi dengan prediktor laten (tidak teramati) disebut regresi laten (Tarpey & Petkova 2010). Salah satu kasus yang dapat diselesaikan dengan model regresi laten yaitu efek plasebo. Efek plasebo adalah efek sugesti yang membuat orang percaya untuk sembuh dengan sejenis obat tertentu. Keunggulan efek plasebo dapat digunakan dalam berbagai hal, terutama dalam pengobatan dan penyembuhan. Salah satu penyakit yang bisa disembuhkan dengan efek plasebo adalah depresi. Depresi merupakan salah satu contoh penyakit yang dapat berubah secara kontinu dan tidak diketahui untuk waktu yang lama. Penelitian yang dilakukan oleh Beecher tentang the powerfull placebo menunjukkan bahwa 32% pasien sembuh karena efek plasebo. Namun, penelitian tersebut tidak memberikan model yang cocok untuk efek plasebo itu sendiri. Oleh karena itu, model regresi laten diusulkan untuk dikaji pada efek plasebo.
Pada penelitian ini, akan dilakukan analisis mengenai model regresi laten pada efek plasebo dengan prediktor laten kontinu menggunakan distribusi beta. Penelitian ini diharapkan dapat memberikan model yang cocok untuk efek plasebo. Metode yang digunakan untuk model regresi laten yaitu algoritma EM (Expectation Maximization). Algoritma EM tergantung pada peubah yang tidak teramati. Algoritma EM dinilai layak digunakan karena terdapat dua tahapan pada setiap iterasi yaitu tahap maximization (M-step) dan tahap expectation (E-step). Algoritma ini mengulang tahapan expectation (E-Step) yang menghitung dugaan kemungkinan dengan memasukkan peubah tersembunyi seolah-olah peubah itu diamati. Tahapan maximization (M-Step) merupakan tahapan yang menghitung dugaan maximum likelihood parameter dengan memaksimalkan dugaan kemungkinan yang diperoleh dari EStep. Nilai parameter yang diperoleh dari MStep digunakan kembali untuk memulai EStep yang selanjutnya. Proses ini akan berulang sehingga dicapai konvergensi nilai likelihood (Dempster et al. 1997). 1.2 Tujuan Tujuan dari karya ilmiah ini adalah untuk menganalisis model regresi laten, menduga parameter dengan metode maximum likelihood dan algoritma EM (Expectation Maximization) serta mengaplikasi model pada efek plasebo dalam pengobatan depresi. 1.3 Sistematika Penulisan Karya ilmiah ini terdiri atas lima bagian. Bagian pertama berupa pendahuluan yang terdiri dari latar belakang, tujuan, dan sistematika penulisan. Bagian kedua merupakan landasan teori yang menyajikan aspek teoritis penulisan karya ilmiah. Bagian ketiga merupakan metode penelitian. Bagian keempat merupakan pembahasan yang membahas model regresi laten beserta aplikasinya untuk efek plasebo. Bagian kelima merupakan simpulan yang diperoleh dari pembahasan karya ilmiah ini.
2
II LANDASAN TEORI 2.1 Analisis Regresi Analisis regresi adalah analisis yang digunakan untuk menganalisis data dan mengambil kesimpulan yang bermakna tentang hubungan kebergantungan yang mungkin ada. Tujuan analisis regresi yaitu untuk mengevaluasi hubungan antara satu peubah dengan satu peubah lainnya atau satu peubah dengan beberapa peubah lainnya. Peubah dapat dibedakan menjadi dua jenis, yaitu peubah prediktor atau peubah bebas, dan peubah respons atau peubah takbebas. Peubah prediktor adalah peubah yang dapat ditentukan atau diatur (misalnya suhu input) atau yang nilainya dapat diamati. Namun, tidak dapat dikendalikan (misalnya kelembaban udara luar). Akibat perubahan yang disengaja atau yang terjadi pada peubah prediktor, suatu pengaruh atau efek dipancarkan ke peubah lain disebut peubah respons. (Draper & Smith 1992) Regresi Linear Sederhana Regresi linear sederhana adalah persamaan regresi yang menggambarkan hubungan antara satu peubah prediktor (x) dan peubah respons (y), hubungan keduanya dinyatakan dalam fungsi linear, sehingga hubungan kedua peubah tersebut dapat dituliskan dalam bentuk persamaan: 𝑦 = 𝛽0 + 𝛽1 𝑥 + 𝜀
(1)
dengan 𝛽0 : parameter regresi, 𝛽1 : parameter regresi, 𝑦 : peubah respons, 𝑥 : peubah prediktor, 𝜀 : galat. (Draper & Smith 1992) Regresi Linear Berganda Regresi linear berganda adalah persamaan regresi yang menggambarkan hubungan antara peubah respons (y) dan banyak peubah prediktor (x). Model regresi linear berganda yang melibatkan p peubah prediktor adalah 𝑦𝑖 = 𝛽0 + 𝛽1 𝑥1𝑖 + 𝛽2 𝑥2𝑖 + ⋯ + 𝛽𝑝 𝑥𝑝𝑖 + 𝜀i.(2)
Bentuk umum regresi berganda adalah p
yi 0 j xij i
(3)
j 1
dengan 𝑦𝑖 : peubah respons, 𝑥𝑖𝑗 : peubah prediktor, 𝛽0 : perpotongan (intercept), 𝛽𝑗 : koefisien parameter dengan j = 1, 2,…, p, 𝜀𝑖 : galat. (Draper & Smith 1992) Pendugaan Koefisien Regresi Linear Berganda Metode kuadrat terkecil adalah suatu metode untuk menghitung koefisien regresi sampel sebagai penduga koefisien regresi populasi (𝜷), sedemikian sehingga jumlah galat kuadrat memiliki nilai terkecil. Secara matematis, model regresi dapat dinyatakan sebagai berikut: 𝒀 = 𝑿𝜷 + 𝛆
(4)
dengan Y : vektor peubah respons berukuran n × m dengan n adalah banyaknya peubah respons yang diamati, X : vektor peubah prediktor berukuran n × p, dengan p adalah banyaknya peubah prediktor, 𝜷 : vektor koefisien berukuran p × m, 𝜺 : vektor galat berukuran n × m yang berdistribusi ℕ(0, 𝝈𝜀 ). (Draper & Smith 1992) Asumsi dasar untuk 𝜀𝑖 atau 𝑦𝑖 , yaitu: 1. 𝐸 𝜺 = 𝟎 atau 𝐸 𝒚 = 𝑿𝜷, 2. 𝐶𝑜𝑣(𝜺) = 𝝈2 𝐈 atau 𝐶𝑜𝑣(𝒚) = 𝝈2 𝐈. Analisis regresi linear berganda menganalisis pengaruh 𝑥1 , 𝑥2 , … , 𝑥𝑝 terhadap y dengan menduga koefisien-koefisien 𝜷0 , 𝜷1 , 𝜷2 , … , 𝜷𝑝 . Analisis ini menggunakan metode jumlah kuadrat terkecil (least sum n
square), yaitu dengan meminimumkan
i 1
2 i
diperoleh nilai dugaan bagi 𝜷. (Rencher & Schaalje 2008)
3
Teorema Estimasi 𝜷 Jika 𝒀 = 𝑿𝜷 + 𝛆, dengan 𝑿𝒏×(𝒑+𝟏) rank 𝑘 + 1 < 𝑛, maka nilai 𝜷 yang n
meminimumkan
i 2 adalah
Jika p menyatakan peluang sukses, maka X merupakan fungsi kerapatan peluang 𝑃𝑋 𝑥 = 𝑝 𝑥 (1 − 𝑝)1−𝑥 , untuk 𝑥 = 0, 1. (Hogg & Craig 1995)
i 1
𝜷 = (𝑿′ 𝑿)−𝟏 (𝑿′ 𝒀).
(5)
Bukti dapat dilihat pada Lampiran 1. (Rencher & Schaalje 2008) 2.2 Peubah Acak dan Distribusi Definisi Peubah Acak Misalkan Ω adalah ruang contoh dari suatu percobaan acak. Fungsi 𝛸 yang terdefinisi pada Ω yang memetakan setiap unsur 𝜔 ∈ Ω ke satu dan hanya satu bilangan real 𝛸 𝜔 = 𝑥 disebut peubah acak. Ruang dari 𝛸 adalah himpunan bagian bilangan real 𝐴 = {𝑥 ∶ 𝑥 = 𝛸 𝜔 , 𝜔 ∈ Ω }. (Hogg et al. 2005) Definisi Peubah Acak Diskret Peubah acak 𝛸 dikatakan diskret jika himpunan semua nilai {𝑥1 , 𝑥2 ,...} merupakan himpunan tercacah. (Grimmett & Stizaker1992) Definisi Peubah Acak Kontinu Peubah acak 𝛸 dikatakan kontinu jika fungsi distribusi komulatifnya 𝐹𝑋 𝑥 adalah fungsi kontinu untuk setiap 𝑥 ∈ 𝑅. (Hogg et al. 2005)
2.4 Family Beta dan Distribusi Beta Definisi Family Beta Misalkan Y suatu peubah acak dengan 𝑓𝑌 (𝑦; 𝑎, 𝑏) merupakan fungsi kepekatan peluang dari peubah acak 𝑌 dengan parameter a dan b, maka 𝑓𝑌 (𝑦; 𝑎, 𝑏) dapat dikatakan sebagai family beta jika dapat dibentuk sebagai berikut: 𝑓𝒀 𝑦 =
1 (𝑦 − 𝑝)𝒂−𝟏 (𝑞 − 𝑦)𝑏−1 (𝑞 − 𝑝)𝒂+𝒃−𝟏 𝐵(𝑎, 𝑏)
(7) dengan 𝑝 ≤ 𝑦 ≤ 𝑞 dan 𝑎 > 0, 𝑏 > 0 merupakan fungsi parameter 𝐵(𝑎, 𝑏). (Johnson et al. 1995) Definisi Distribusi Beta Suatu peubah acak 𝑋 dikatakan mempunyai distribusi beta dengan parameter 𝛼 > 0 dan 𝛽 > 0 jika fungsi kepekatannya diberikan oleh f X ( x)
1 x 1 (1 x) 1 B( , )
(8)
1
dengan B( , ) x 1 (1 x) 1 dx dan 0
Definisi Fungsi Distribusi Jika 𝑋 suatu peubah acak, distribusinya didefinisikan sebagai 𝐹𝑋 𝑥 = 𝑃[𝑋 ≤ 𝑥]
0 < 𝑥 < 1. Rataan dan ragam : fungsi
(6)
untuk setiap 𝑥 ∈ (−∞, +∞). (Ghahramani 2000) 2.3 Distribusi Bernoulli Suatu percobaan acak yang menghasilkan dua kemungkinan (sukses dan gagal) disebut percobaan Bernoulli. Peubah acak X disebut mempunyai sebaran Bernoulli, jika X merupakan peubah acak pada percobaan
1, jika sukses, Bernoulli dengan X 0, jika gagal.
2
(9)
,
(10) . ( ) ( 1) (Hogg & Craig 1995) 2
2.5 Matriks Definisi Matriks Identitas Matriks identitas adalah matriks 𝐼 = 𝛿𝑖𝑗 yang berorde 𝑛 × 𝑛, dengan 𝛿𝑖𝑗 =
1, 0,
𝑖 = 𝑗, 𝑖 ≠ 𝑗.
(11) (Leon 1998)
4
Definisi Matriks Simetris Suatu matriks A berorde 𝑛 × 𝑛 disebut simetris jika 𝐴′ = 𝐴. (Leon 1998) Definisi Invers Suatu matriks 𝐴 yang berorde 𝑛 × 𝑛 dikatakan taksingular jika terdapat matriks 𝐵 sehingga 𝐴𝐵 = 𝐵𝐴 = 𝐼. Matriks 𝐵 dikatakan invers multiplikatif dari matriks 𝐴. Invers multiplikatif dari matriks taksingular 𝐴 secara sederhana disebut juga sebagai invers dari matriks 𝐴 dan dinotasikan dengan 𝐴−1 . (Leon 1998) Definisi Transpos Transpos dari suatu matriks 𝐴 = 𝑎𝑖𝑗 berorde 𝑚 × 𝑛 adalah matriks 𝐵 = 𝑏𝑗𝑖 berorde 𝑛 × 𝑚 yang didefinisikan oleh 𝑏𝑗𝑖 = 𝑎𝑖𝑗 untuk setiap 𝑖 dan 𝑗. Transpos dari 𝐴 dinotasikan oleh 𝐴′ . (Leon 1998) Definisi Trace Trace adalah fungsi yang didefinisikan hanya pada matriks persegi. Misalkan A adalah matriks 𝑛 × 𝑛, trace dari matriks A adalah jumlah dari elemen-elemen diagonal dari matriks A, 𝑛
𝑡𝑟 𝐴 =
𝑎𝑖𝑖 . 𝑖=1
(12)
(Magnus & Neudecker 1999) Teorema Trace Misalkan 𝜆 adalah skalar dan A dan B adalah matriks yang masing-masing memiliki orde 𝑛 × 𝑛, maka:
1. 2. 3. 4.
tr(𝐴 + 𝐵) = tr 𝐴 + tr 𝐵, tr(𝜆𝐴) = 𝜆 tr(𝐴), tr(𝐴′) = tr(𝐴), tr(𝐴𝐵) = tr(𝐵𝐴).
Bukti dapat dilihat pada Lampiran 4. (Magnus & Neudecker 1999) 2.6 Metode Maximum Likelihood Misalkan 𝑋1 , 𝑋2 , … , 𝑋𝑛 adalah peubah acak i.i.d dengan fungsi kepekatan peluang 𝑓(𝑥; 𝜃), dengan 𝜃 diasumsikan skalar dan tidak diketahui, maka prosedur fungsi likelihood dapat dituliskan sebagai berikut: 𝑛
𝐿 𝑥, 𝜃 =
𝑓(𝑥𝑖 , 𝜃) , 𝜃 ∈ Ω 𝑖=1
(13)
dengan 𝑋 = (𝑋1 , 𝑋2 , … , 𝑋𝑛 ). Fungsi loglikelihood dari 𝐿 𝑥, 𝜃 , dapat dinotasikan dengan: 𝑙 𝑥, 𝜃 = log 𝐿 𝑥, 𝜃 𝑛
=
log 𝑓 𝑥𝑖 , 𝜃 , 𝜃 ∈ Ω . 𝑖=1
(14) Pendugaan parameter dengan metode maximum likelihood estimation dapat 𝜕𝑙 𝑥,𝜃 diperoleh dari = 0. 𝜕𝜃 (Hogg & Craig 1995) 2.7 Algoritma K-Means Algoritma K-Means merupakan salah satu metode pengelompokan data tidak berhirarki yang berusaha mempartisi data yang ada ke dalam satu atau lebih gerombol. (Hair et al. 1995)
III METODE PENELITIAN Karya ilmiah ini disusun melalui studi literatur mengenai regresi laten, sumber data pasien depresi, analisis dan pemrograman dari data sekunder menggunakan software R 2.13.1. 3.1 Studi Literatur Studi literatur meliputi pencarian berbagai informasi yang berhubungan dengan topik yang dibahas yaitu regresi laten.
Langkah-langkah penelitian: 1. Menelusuri model regresi laten yang sesuai, 2. Mengestimasi parameter model regresi ketika x berupa laten menggunakan Kmeans, maximum likelihood dan algoritma EM, 3. Melakukan aplikasi pada efek plasebo.
5
3.2 Sumber Data Data yang digunakan untuk penelitian ini adalah data sekunder mengenai pasien depresi rawat jalan yang bersumber dari Tarpey dan Petkova (2010) dengan komunikasi via email. Objek penelitian adalah pasien laki-laki dan perempuan berusia 18 - 65 tahun. Jumlah responden depresi dalam penelitian adalah 393 orang. Tarpey dan Petkova melihat perubahan tingkat gejala depresi menggunakan skala Hamilton Depression Rating (HAM-D). 3.3 Analisis dan Pemrograman Tahapan analisis dan pemrograman yang dilakukan dalam penelitian ini sebagai berikut:
1. Membangun model regresi laten Tahap pertama ini akan dibangun model regresi dengan inisialisasi setiap parameter, kemudian disimpan dalam bentuk file emnorm.r. 2. Melakukan pengelompokan data dengan algoritma K-means, yaitu dengan menentukan inisialisasi nilai awal k, komponen means dan perkiraan matriks kovarian. 3. Melakukan pendugaan algoritma EM, yaitu dengan menentukan inisialisasi nilai sigma, perkiraan nilai parameter beta, mengerjakan tahap E-step dan M-step.
IV PEMBAHASAN 4.1 Efek Plasebo Plasebo adalah istilah medis untuk sejenis obat tanpa bahan kimia yang kadang hanya berisi gula atau cairan garam. Efek plasebo adalah efek sugesti yang membuat orang percaya untuk sembuh dengan sejenis obat tertentu. Efek plasebo digunakan dalam pengobatan dan penyembuhan. Efek ini lebih menekankan pada faktor psikologis dan keyakinan untuk sembuh. 4.2 Model Regresi Laten Analisis multivariat merupakan salah satu jenis analisis statistik yang digunakan untuk menganalisis data dengan data yang digunakan berupa banyak peubah prediktor dan banyak peubah respons. Bentuk hubungan analisis regresi linear ganda yaitu beberapa peubah prediktor terhadap satu peubah respons. Model regresi laten pada efek plasebo: 𝒀 = 𝑿𝜷 + 𝜺
(15)
dengan 𝒀 : peubah respons p-variat, 𝑿 : prediktor laten (tidak teramati) berdimensi r, 𝜷 : parameter, 𝜺 : galat (𝜺 ~ ℕ(0,𝚿)). Penjelasan model regresi linear multivariat dapat dilihat pada Lampiran 2. Jika populasi terdiri dari dua sub populasi laten (orang yang mengalami efek plasebo dan orang yang tidak
mengalami efek plasebo), maka prediktor laten pada persamaan (15) memiliki distribusi Bernoulli. Misalkan p = P(x = 1) dan asumsikan bahwa 𝜀 menyebar distribusi ℕ(𝟎, 𝚿) dengan 𝚿 adalah matriks kovarian definit positif. Kepadatan marjinal 𝑓(𝒚) pada y dinamakan model campuran terbatas: 𝑓 𝒚 =
𝑓(𝑥, 𝒚) 𝑥
= 𝑃 𝑥 = 0 𝑓 𝒚 𝑥 = 0 +𝑃 𝑥 = 1 𝑓(𝒚|𝑥 = 1) = 1 − 𝑝 ℕ 𝒚; 𝝁1 , 𝚿 + 𝑝ℕ(𝒚; 𝝁2 , 𝚿) (16) dengan 𝝁1 = 𝛽0 , 𝝁2 = 𝛽0 + 𝛽1 dan ℕ(𝒚; 𝝁, 𝚿) merupakan hasil fungsi kepadatan normal berganda dengan rataan vektor 𝝁 dan matriks kovarians 𝚿. Distribusi Bernoulli memiliki probabilitas 0 dan 1. Cara alami untuk menggeneralisasi model campuran terbatas yaitu dengan mengganti prediktor Bernoulli 0 − 1 oleh distribusi kontinu pada interval (0, 1) yang kepadatannya berbentuk U. Pilihan alami untuk distribusi pada x adalah distribusi beta dengan kepadatan 𝑔(𝑥; 𝑎, 𝑏) yang didefinisikan dalam parameter a dan b. 𝑔 𝑥; 𝑎, 𝑏 =
Γ 𝑎+𝑏 Γ 𝑎 Γ 𝑏
𝑥 𝑎−1 1 − 𝑥
𝑏−1
(17) dengan 0 < 𝑥 < 1. Family beta menghasilkan berbagai bentuk kepadatan salah satunya kepadatan berbentuk U. Kepadatan ini menyediakan generalisasi kontinu dari distribusi Bernoulli diskret.
6
Kepadatan bersama untuk x dan y dalam persamaan (15) adalah 𝑓 𝑥, 𝒚 = 𝑓 𝒚 𝑥; 𝛽0 , 𝛽1 , 𝚿 𝑔 𝑥; 𝑎, 𝑏 = ℕ(𝒚; 𝛽0 + 𝛽1 𝑥, 𝚿)𝑔 𝑥; 𝑎, 𝑏
(18)
dengan 𝑔 𝑥; 𝑎, 𝑏 merupakan distribusi beta yang diberikan dalam persamaan (17). Kepadatan marjinal hasilnya adalah 1
𝑓 𝒚 = 0
1 2𝜋𝜎
exp{−(𝒚 − 𝛽0 − 𝛽1 𝑥)′𝚿−1
𝒚 − 𝛽0 − 𝛽1 𝑥 /2}𝑔(𝑥; 𝑎, 𝑏)𝑑𝑥. (19)
4.3 Metode Pendugaan Parameter Ada beberapa metode yang digunakan untuk menduga parameter, antara lain metode momen, metode Bayes, metode maximum likelihood, dan algoritma EM (ExpectationMaximization). Pendugaan parameter yang digunakan pada penelitian ini adalah metode maximum likelihood dan algoritma EM (Expectation-Maximization). Metode Maximum Likelihood Metode maximum likelihood adalah suatu metode yang baik untuk memperoleh sebuah parameter tunggal. Misalkan 𝑋1 , 𝑋2 , … , 𝑋𝑛 masing-masing peubah acak saling bebas dengan sebaran yang memiliki fungsi kepekatan peluang 𝑓 𝑥; 𝜃 , dengan 0 ≤ 𝜃 ≤ 1, 𝜃 𝜖 Ω dan Ω adalah ruang contoh. Fungsi kepekatan peluang bersama dari 𝑋1 , 𝑋2 , … , 𝑋𝑛 adalah 𝐿 𝜃 𝑥1 , 𝑥2 , … , 𝑥𝑛 = 𝑓 𝑥1 𝜃 𝑓 𝑥2 𝜃 … 𝑓 𝑥𝑛 𝜃 yang disebut juga sebagai fungsi likelihood. Fungsi sederhana dari x1, x2, … ,xn yaitu (x1, x2, … , xn), sehingga 𝜃 = u(x1, x2, … , xn) membuat fungsi kemungkinan L maksimum untuk semua 𝜃 𝜖 Ω. Statistik u(x1, x2, … , xn) disebut penduga maximum likelihood dari 𝜃 yang dinotasikan dengan 𝜃 = u(X1, X2, … , Xn). Untuk menduga parameter dengan menggunakan metode maximum likelihood tidak bisa secara langsung karena datanya tidak teramati, untuk itu dapat digunakan algoritma EM. Algoritma Expectation-Maximization Algoritma Expectation-Maximization adalah suatu algoritma yang sangat handal untuk pendugaan parameter dengan menggunakan metode maximum likelihood dari fungsi likelihood pada data yang tidak teramati (McLachlan & Krishnan 2008). Proses algoritma EM dilakukan dengan dua tahap, yaitu:
E-Step (Tahapan Expectation) Merupakan tahapan untuk menghitung ekspektasi bersyarat dari fungsi loglikelihood dengan prediktor laten. Misalkan 𝜷 adalah suatu nilai awal, maka E-Step didefinisikan 𝐸(𝑙𝑛𝐿(𝜷; 𝚿)|𝒀)
(20)
dalam aplikasi pada efek placebo, didefinisikan sebagai matriks koefisien: 𝜷=
𝜷0 ′ 𝜷1 ′
𝜷
(21)
dengan 𝜷 berdimensi 2 × p, masing-masing kolom p dari 𝜷 merupakan intercept dan slope koefisien regresi untuk setiap peubah respons p. X merupakan matriks yang kolom pertama terdiri dari 1 dan kolom kedua terdiri dari prediktor laten 𝑥𝒊, 𝑖 = 1, 2, …, n, Y merupakan matriks berdimensi n × p, dan 𝚿 merupakan matriks kovarian definit positif. Model regresi laten dapat ditulis sebagai 𝒀 = 𝑿𝜷 + 𝜺. Likelihood untuk model dinyatakan sebagai berikut : 𝑛𝑝
(22) regresi
laten
𝑛
𝐿 𝜷, 𝚿 = 2𝜋 − 2 𝚿 −2 exp{−𝑡𝑟[𝚿−1 𝒀 − 𝑿𝜷 ′ (𝒀 − 𝑿𝜷)]/2}. (23) Misalkan 𝑿 = 𝐸[𝑿|𝒀] dan 𝑿′ 𝑿 = 𝐸[𝑿′𝑿|𝒀], maka nilai harapan dari 𝑡𝑟[𝚿−1 𝒀 − 𝑿𝜷 ′ (𝒀 − 𝑿𝜷) bersyarat 𝒀 adalah 𝐸 𝑡𝑟 𝚿−1 𝒀 − 𝑿𝜷 ′ 𝒀 − 𝑿𝜷 𝒀 = 𝑡𝑟[𝚿−1 {𝒀′ 𝒀 − 𝒀′ 𝑿𝜷 − 𝜷′𝑿′𝒀 + 𝜷′ (𝑿′𝑿)𝜷}]. (24) Pembuktian persamaan (24) dapat dilihat pada Lampiran 3. 𝑿′ 𝑿 dan 𝚿 definit positif, sehingga pada E-Step diperoleh 𝑛𝑝 𝑛 ln 2𝜋 − ln 𝚿 2 2 −0.5𝑡𝑟[𝚿−1 {𝒀′ 𝒀 − 𝒀′ 𝑿𝜷 − 𝜷′𝑿′𝒀 +𝜷′ (𝑿′𝑿)𝜷}]. (25)
𝐸(𝑙𝑛𝐿(𝜷; 𝚿)|𝒀) = −
Misalkan 𝜷 = (𝑿′ 𝑿)−1 𝑿′𝒀 maka persamaan (24) dapat dinyatakan kembali sebagai: 𝑡𝑟[𝚿−1 {(𝒀′ 𝒀 − 𝜷′ 𝑿′ 𝑿 𝜷) + 𝜷 − 𝜷 𝑿′ 𝑿 (𝜷 − 𝜷)}].
′
(26)
7
Latent Beta Distribution for Simulation
= 𝑡𝑟 𝑿′ 𝑿 𝜷 − 𝜷 𝚿−1 𝜷 − 𝜷
8 6
𝑿′ 𝑿 𝜷 − 𝜷
Density
′
′
4
𝑡𝑟 𝚿 −1 𝜷 − 𝜷
10
Pembuktian persamaan (26) dapat dilihat pada Lampiran 4. Bagian dari trace pada persamaan (26) yang melibatkan parameter 𝜷 dapat dinyatakan sebagai berikut:
′
Hasil E-Step sebagai:
dapat
dinyatakan
(27) kembali
𝑛 𝑛𝑝 𝐸(𝑙𝑛𝐿(𝜷; 𝚿)|𝒀) = − ln 2𝜋 − ln 𝚿 2 2 −0.5𝑡𝑟[𝚿−1 {𝒀′ 𝒀 − 𝜷′ 𝑿′ 𝑿 𝜷 − 𝜷′ 𝑿′ 𝑿 𝜷 +𝜷′ (𝑿′𝑿)𝜷}]. (28) M-Step (Tahapan Maximization) Merupakan tahapan untuk mendapatkan parameter baru 𝜷 dengan memaksimumkan 𝐸(𝑙𝑛𝐿(𝜷; 𝚿)|𝒀), yang dinyatakan sebagai berikut: 𝜕(𝐸(𝑙𝑛𝐿(𝜷; 𝚿)|𝒀)) 𝜕𝜷 sehingga diperoleh 𝜷 = 𝜷 = 𝑿′ 𝑿
(29) −1
𝑿′ 𝒀.
Proses E-Step dan M-Step ini dilakukan terus secara iteratif sampai diperoleh suatu nilai dugaan parameter 𝜷 yang konvergen. Langkah-langkah mencari nilai dugaan parameter 𝜷 menggunakan algoritma EM dapat dilihat pada Lampiran 5.
2
= 𝑡𝑟[(𝑿′ 𝑿)1/2 (𝜷 − 𝜷)𝚿−1 𝜷 − 𝜷 (𝑿′ 𝑿)1/2 ].
0.0
0.2
0.4
0.6
0.8
1.0
xx
Gambar 1 Plot kepadatan beta dengan a = 0.5 dan b = 0.3 untuk peubah laten x. Kurva pada Gambar 1 merupakan kurva distribusi beta dengan x peubah laten. Kurva tersebut menghasilkan kerapatan berbentuk U. Hasilnya mendukung bahwa tidak ada dua kelas berbeda dari subjek (orang-orang yang mengalami efek plasebo dan orang-orang yang tidak mengalaminya). Kode R untuk Gambar 1 dapat dilihat pada lampiran 7. Jumlah responden depresi dalam Tarpey dan Petkova sebanyak 393 orang dalam satu minngu. Salah satu cara yang digunakan dalam menduga sebaran pasien depresi selama satu minggu adalah dengan melihat kesesuaian histogram data. Histogram data pasien depresi rawat jalan selama satu minggu dapat dilihat pada Gambar 2.
4.4 Aplikasi pada Efek Plasebo Algoritma EM yang dijelaskan terdahulu akan diuji menggunakan berbagai pengaturan parameter untuk model regresi laten. Untuk setiap pengaturan parameter, 50 himpunan data diberikan dengan masing-masing ukuran sampel n = 100. Sebagai ilustrasi, peubah laten x menyebar distribusi beta dengan parameter a = 0.5 dan b = 0.3 menghasilkan kerapatan berbentuk U. Galat menyebar distribusi normal dengan rataan 𝜇 = 0 dan standar deviasi 𝜎 = 0.5. Gambar 2 Histogram data pasien depresi rawat jalan.
8
Histogram pada Gambar 2 menunjukkan bahwa data pasien depresi berdistribusi beta. Pasien depresi yang mendapatkan efek plasebo ternyata mengalami perubahan positif. Efek plasebo telah diketahui ada di penelitian depresi dan efek ini cenderung untuk meningkatkan suasana hati. Tampaknya masuk akal atribut peubah laten untuk efek plasebo dalam penelitian ini. Karena efek plasebo adalah laten, model dari gejala depresi pada minggu pertama merupakan fungsi kekuatan efek plasebo x.
Jika x dalam persamaan (15) adalah Bernoulli, maka model menjelaskan bahwa terdapat dua jenis subjek. Dua jenis subjek tersebut, yaitu orang-orang yang mengalami efek plasebo dan orang-orang yang tidak mengalaminya. Jika x kontinu, maka model menjelaskan bahwa masing-masing subjek mengalami efek plasebo berderajat variasi. Perubahan tingkat gejala depresi dari awal sampai minggu pertama diukur pada skala Hamilton Depression Rating (HAM-D).
Non-parametrik --- Laten --- Campuran
Gambar 3 Histogram penurunan depresi. Perubahan positif dari awal sampai minggu pertama menunjukkan penurunan depresi seperti terlihat pada Gambar 3. Kurva pada Gambar 3 merupakan suatu perhitungan kerapatan non-parametik dari y yang menunjukkan distribusi miring (condong) kanan. Perubahan positif menunjukkan berbagai derajat peningkatan antara subjek yang relatif baik. Kurva solid merupakan perkiraan kepadatan non-parametrik, kurva putus-putus merah merupakan perkiraan kepadatan regresi laten, dan kurva hijau putusputus merupakan perkiraan kepadatan campuran dua komponen terbatas.
Estimasi model regresi laten adalah 𝑦 = −0.279 + 26.604 𝑥 yang ditunjukkan pada Gambar 4. Interpretasinya adalah semakin besar nilai efektivitas plasebo maka akan memberikan perubahan suasana hati yang semakin tinggi. Pada awalnya (minggu pertama) efek plasebo berubah secara kontinu, mulai dari yang sangat lemah sampai cukup kuat sehingga distribusi efek plasebo condong ke arah lebih baik. Efek plasebo masih bervariasi tetapi kita mulai melihat segmentasi ke dalam dua kelompok.
Latent Regression Quantiles
25
y Hati Perubahan Suasana
20
15
10
𝑦 = −0.279 + 26.604 𝑥
5
0 0.0
0.2
0.4 0.6 x (Tingkat Efektivitas Plasebo)
0.8
1.0
x
Gambar 4 Grafik estimasi model regresi laten.
Gambar 5 Plot kuantil-kuantil model regresi laten untuk data depresi.
9
Plot kuantil dari histogram pada minggu pertama dapat dilihat pada Gambar 5. Berdasarkan Gambar 5, pola pancaran titiktitik untuk model regresi laten hampir membentuk garis lurus tetapi pada ekor kanan menyimpang dari rentang data. Hal ini dikarenakan terdapat dua kelompok, yaitu kelompok yang tidak ada pengaruh efek plasebo dan kelompok yang kuat pengaruh plasebo. Berdasarkan ide dasar uji SharpioWilks untuk normalitas, 𝑅2 yang diperoleh sebesar 0.95. Koefisien 𝑅2 ini menunjukkan
bahwa pendekatan model regresi laten cocok untuk gugus data pengamatan efek placebo pada studi depresi. Pendugaan parameter model dilakukan dengan menggunakan Software-R. Kode R sebagai pendefinisian awal dapat dilihat pada Lampiran 8 dan kode R untuk algoritma Kmeans, maximum likelihood dan algoritma EM dapat dilihat pada Lampiran 9. Kode R untuk plot kuantil-kuantil model regresi laten untuk data depresi dapat dilihat pada Lampiran 10.
V SIMPULAN DAN SARAN 5.1 Simpulan Simpulan yang dapat diambil setelah melakukan penelitian dan pengumpulan data tentang model matematika, yaitu semakin besar nilai efektivitas plasebo maka akan memberikan perubahan suasana hati yang semakin tinggi. Persamaan regresi laten yang menggambarkan efek plasebo pada studi depresi adalah 𝑦 = -0.279 + 26.604 x. Parameter model regresi laten tersebut diduga dengan menggunakan metode maximum
likelihood dan algoritma EM (ExpectationMaximization). 5.2 Saran 1. Penelitian lebih lanjut diharapkan mengkaji dengan distribusi lain yang memungkinkan untuk prediktor laten dengan faktor-faktor lain yang memengaruhi gejala depresi pada efek plasebo. 2. Perlu diketahui substansi dari fenomena masalah yang dikaji untuk memilih model yang dianggap cocok, sehingga pemilihan model akan sesuai dengan fakta yang dijumpai.
DAFTAR PUSTAKA Dempster AP, Laird NM, Rubin D. 1997. Maximum likelihood from incomplete data via the EM algorithm. Journal of the Royal Statistical Society. B 39, 1-38. Draper N, Smith H. 1992. Analisis Regresi Terapan. Ed ke-2. Jakarta: PT Gramedia Pustaka Utama. Ghahramani S. 2005. Fundamental of Probability and Random Processes. Ed ke2. New Jersey: Prentice Hall. Grimmett GR, Stirzaker DR. 1992. Probability and Random Processes. Ed ke2. Oxford: Clarendon Press. Hair JF et al. 1995. Multivariate Data Analysis. Ed ke-4. New Jersey: Prentice Hall. Hogg RV, Craig AT. 1995. Introduction to Mathematical Statistics. Ed ke-5. New Jersey: Prentice Hall. Hogg RV, Craig AT, McKean JW. 2005. Introduction to Mathematical Statistics. Ed ke-6. New Jersey: Prentice Hall.
Johnson NL, Kotz S, Balakrishnan N. 1995. Continuous Univariate Distributions. New York: J Wiley. Leon SJ. 1998. Linear Algebra with Applications. Ed ke-5. New Jersey: Prentice Hall. Magnus JR, Neudecker. 1999. Matrix Differential Calculus with applications in statistics and Econometrics. Ed revisi. New York: J Wiley. McLachlan GJ, Krishnan T. 2008. The EM Algorithm and Extension. Ed ke-2. New Jersey: J Wiley. Rencher AC, Schaalje GB. 2008. Linear Models in Statistics. Ed ke-2. New Jersey: J Wiley. Tarpey T, Petkova E. 2010. Latent regression analysis. Statistical Modelling. 10(2): 133158.
10
LAMPIRAN
11
Lampiran 1 Pembuktian Teorema Estimasi 𝜷 Bentuk umum regresi linear berganda:
𝑝
𝑦𝑖 = 𝛽0 +
𝛽𝑗 𝑥𝑖𝑗 + 𝜀𝑖 𝑗 =1
Persamaan diobservasi sebanyak n pengamatan sehingga: y1 0 1 x11 2 x12 ... p x1 p 1 x x ... x 1 21 2 22 p 2p 2 y2 0 yn 0 1 xn1 2 xn 2 ... p xnp n 1 x11 x12 ... x1 p 0 1 1 x21 x22 ... x2 p 1 2 1 x x ... x p n np n1 n 2
Metode Kuadrat Terkecil (Least Square) n
Meminimumkan
i 1
n
i 1
j dengan j 1, 2,..., p , ideal n
i 1
n
2 i
2
( yi yi ) dengan kendala parameternya, yaitu:
2 i
n
i 1
2 i
0 0
^
( yi yi )2 i 1 n
p
i 1
j 1
[ yi ( ˆ0 j xij )]2 n
p
i 1
j 1
p
[ yi 2 2 yi ( ˆ0 j xij ) ( ˆ0 j xij ) 2 ] j 1
n
p
p
p
i 1
j 1
j 1
j 1
[ yi 2 2 yi ˆ0 2 yi j xij ˆ0 2 2ˆ0 j xij j 2 xij 2 ] n
( i 2 ) i 1
j
p
n
0 (2 xij yi 2ˆ0 xij 2 j xij 2 ) 0 i 1 j 1
p
n
2 ( xij yi ˆ0 xij j xij 2 ) 0 i 1 j 1
n
p
( xij yi ˆ0 xij j xij 2 ) 0 i 1 j 1
Dalam bentuk matriks: 𝛆′𝛆 = (𝒀 − 𝑿𝜷)′ (𝒀 − 𝑿𝜷) = 𝐘 ′ 𝐘 − 𝜷′ 𝑿′ 𝒀 − 𝒀′ 𝑿𝜷 + 𝜷′ 𝑿′ 𝑿𝜷 = 𝒀′ 𝒀 − 𝟐𝜷′ 𝑿′ 𝒀 + 𝜷′ 𝑿′ 𝑿𝜷
12
Lampiran 1 (Lanjutan)
𝝏(𝛆′𝛆) =𝟎 𝝏𝜷 ′ ′ ′ ′ ′ 𝝏(𝒀 𝒀 − 𝟐𝜷 𝑿 𝒀 + 𝜷 𝑿 𝑿𝜷) =𝟎 ↔ 𝝏𝜷 ↔ −𝟐𝑿′ 𝒀 + 𝟐𝑿′ 𝑿𝜷 = 𝟎 ↔ −𝑿′(𝒑 × 𝒏) 𝒀(𝒏 × 𝟏) + 𝑿′(𝒑 × 𝒏) 𝑿𝒏 × 𝒑 𝜷(𝒒 × 𝟏) = 𝟎 ↔ 𝑿′ 𝒀 = 𝑿′𝑿𝜷 ↔ 𝑿′ 𝑿 −𝟏 𝑿′ 𝑿 𝜷 = 𝑿′ 𝑿 −𝟏 𝑿′ 𝒀 ↔ 𝑰𝜷 = 𝑿′ 𝑿 −𝟏 𝑿′ 𝒀 ↔ 𝜷 = 𝑿′ 𝑿 −𝟏 𝑿′ 𝒀
Sehingga diperoleh 𝜷 = 𝑿′ 𝑿 Bukti lengkap. ∎
−𝟏
𝑿′ 𝒀.
13
Lampiran 2 Penjelasan Model Regresi Linear Multivariat Misalkan untuk peubah respons sebanyak 2 atau terdapat Y1 dan Y2 serta 3 peubah prediktor, maka 𝑦11 𝑦21 ⋮ 𝑦𝑛1
1 𝑥11 𝑥12 𝑦12 1 𝑥21 𝑥22 𝑦22 = ⋮ ⋮ ⋮ ⋮ 𝑦𝑛2 1 𝑥𝑛1 𝑥𝑛2
𝑥13 𝑥23 ⋮ 𝑥𝑛3
𝛽01 𝛽11 𝛽21 𝛽𝑛1
𝜀11 𝛽02 𝜀21 𝛽12 + ⋮ 𝛽22 𝜀𝑛1 𝛽𝑛2
𝜀12 𝜀22 ⋮ . 𝜀𝑛2
Jika terdapat p peubah respons Y dan r peubah prediktor x, terdapat sejumlah persamaan model regresi: 𝑦11 ⋮ 𝑦𝑛1
⋯ ⋱ ⋯
1 𝑥11 𝑥12 𝑦1𝑝 𝑥 ⋮ = 1 𝑥21 22 ⋮ ⋮ ⋮ 𝑦𝑛𝑝 1 𝑥𝑛1 𝑥𝑛2
𝑥1𝑟 𝑥23 ⋮ 𝑥𝑛𝑟
𝛽01 𝛽11 ⋮ 𝛽𝑟1
𝛽02 𝛽12 ⋮ 𝛽𝑟2
𝜀11 … 𝛽0𝑝 𝜀21 … 𝛽1𝑝 + ⋮ ⋱ ⋮ 𝜀𝑛1 … 𝛽𝑟𝑝
dengan 𝜺 = [𝜀1 , 𝜀2 , … , 𝜀𝑝 ], 𝐸 𝜺 = 0, Var(𝜺) = 𝚺. Model regresi linear multivariat adalah 𝒀(𝑛×𝑝) = 𝑿(𝑛×(𝑟+1)) 𝜷((𝑟+1)×𝑝) + 𝜺 𝑛×𝑝 . Lengkap. ∎
𝜀1𝑝 𝜀2𝑝 ⋮ 𝜀𝑛𝑝
14
Lampiran 3 Pembuktian Persamaan (24) Diketahui 𝐸[𝑡𝑟[𝚿−1 𝒀 − 𝑿𝜷 ′ (𝒀 − 𝑿𝜷)]|𝒀] Akan dibuktikan Persamaan (22), yaitu 𝑡𝑟[𝚿−1 {𝒀′ 𝒀 − 𝒀′ 𝑿𝜷 − 𝜷′𝑿′𝒀 + 𝜷′ (𝑿′𝑿)𝜷}] dengan 𝐸 𝑿𝒀 =𝑿 𝐸 𝑿′ 𝑿 𝒀 = 𝑿′ 𝑿 Bukti: 𝐸[𝑡𝑟[𝚿−1 𝒀 − 𝑿𝜷 ′ (𝒀 − 𝑿𝜷)]|𝒀] ↔ 𝐸[𝑡𝑟[𝚿−1 (𝒀′ − 𝜷′ 𝑿′)(𝒀 − 𝑿𝜷)]|𝒀] ↔ 𝐸[𝑡𝑟 𝚿 −1 𝒀′ 𝒀 − 𝒀′ 𝑿𝜷 − 𝜷′ 𝑿′ 𝒀 + 𝜷′ 𝑿′ 𝑿𝜷 𝒀 ] ↔ 𝐸[𝑡𝑟 𝚿 −1 𝒀′ 𝒀 − 𝒀′ 𝑿𝜷 − 𝜷′ 𝑿′ 𝒀 + 𝜷′ 𝑿′ 𝑿 𝜷 𝒀 ] ↔ 𝑡𝑟[𝚿−1 {𝒀′ 𝒀 − 𝒀′ 𝑿𝜷 − 𝜷′𝑿′𝒀 + 𝜷′ (𝑿′𝑿)𝜷}] Bukti lengkap. ∎
15
Lampiran 4 Pembuktian Persamaan (26) Diketahui Persamaan (22), yaitu 𝑡𝑟[𝚿−1 {𝒀′ 𝒀 − 𝒀′ 𝑿𝜷 − 𝜷′𝑿′𝒀 + 𝜷′ (𝑿′𝑿)𝜷}] Akan dibuktikan bahwa Persamaan (24) ekivalen dengan Persamaan (22), yaitu 𝑡𝑟[𝚿−1 {𝒀′ 𝒀 − 𝜷′ 𝑿′ 𝑿 𝜷 + (𝜷 − 𝜷)′ 𝑿′ 𝑿 (𝜷 − 𝜷)}] 𝑿 merupakan matriks yang kolom pertama terdiri dari 1 dan kolom kedua terdiri dari prediktor laten 𝑥𝒊, i = 1, 2, … , n sehingga 𝑿′ 𝑿 merupakan matriks ukuran 2 × 2 yang simetri, maka: 𝑿′ 𝑿 ′ = 𝑿′ 𝑿 dengan 𝜷 = (𝑿′𝑿)−1 𝑿′𝒀 ↔ 𝑿′ 𝑿 𝜷 = 𝑿′ 𝑿 (𝑿′𝑿)−1 𝑿′𝒀 ↔ 𝑿′ 𝑿 𝜷 = 𝑿′𝒀 ↔ ( 𝑿′ 𝑿 𝜷)′ = (𝑿′ 𝒀)′ ↔ 𝜷′ 𝑿′ 𝑿 ′ = 𝒀′ 𝑿 ↔ 𝜷′ 𝑿′ 𝑿 = 𝒀′ 𝑿 Bukti: 𝑡𝑟[𝚿−1 {𝒀′ 𝒀 − 𝒀′ 𝑿𝜷 − 𝜷′𝑿′𝒀 + 𝜷′ (𝑿′𝑿)𝜷}] ↔ 𝑡𝑟[𝚿−1 {𝒀′ 𝒀 − 𝜷′ 𝑿′ 𝑿 𝜷 − 𝜷′𝑿′𝒀 + 𝜷′ (𝑿′𝑿)𝜷}] ↔ 𝑡𝑟[𝚿−1 {𝒀′ 𝒀 − 𝜷′ 𝑿′ 𝑿 𝜷 − 𝜷′ 𝑿′ 𝑿 𝜷 + 𝜷′ (𝑿′𝑿)𝜷}] ↔ 𝑡𝑟[𝚿−1 {𝒀′ 𝒀 − 𝜷′ 𝑿′ 𝑿 𝜷 − 𝜷′ 𝑿′ 𝑿 𝜷 + 𝜷′ (𝑿′ 𝑿)𝜷 − 𝜷′ 𝑿′ 𝒀 + 𝜷′ 𝑿′ 𝒀}] ↔ 𝑡𝑟[𝚿−1 {𝒀′ 𝒀 − 𝜷′ 𝑿′ 𝑿 𝜷 − 𝜷′ 𝑿′ 𝑿 𝜷 + 𝜷′ (𝑿′ 𝑿)𝜷 − 𝜷′ 𝑿′ 𝑿 𝜷 + 𝜷′ 𝑿′ 𝑿 𝜷}] ↔ 𝑡𝑟[𝚿−1 {𝒀′ 𝒀 − 𝜷′ 𝑿′ 𝑿 𝜷 + 𝜷′ 𝑿′ 𝑿 𝜷 − 𝜷′ 𝑿′ 𝑿 𝜷 − 𝜷′ 𝑿′ 𝑿 𝜷 + 𝜷′ 𝑿′ 𝑿 𝜷}] ↔ 𝑡𝑟[𝚿−1 {𝒀′ 𝒀 − 𝜷′ 𝑿′ 𝑿 𝜷 + (𝜷′ 𝑿′ 𝑿 − 𝜷′ 𝑿′ 𝑿 )(𝜷 − 𝜷)}] ↔ 𝑡𝑟[𝚿−1 {𝒀′ 𝒀 − 𝜷′ 𝑿′ 𝑿 𝜷 + (𝜷 − 𝜷)′ 𝑿′ 𝑿 (𝜷 − 𝜷)}] Bukti lengkap. ∎
16
Lampiran 5 Pembuktian Teorema Trace Diketahui trace:
𝑛
𝑡𝑟 𝐴 =
𝑎𝑖𝑖 𝑖=1
Misalkan 𝜆 adalah skalar dan A dan B adalah matriks yang masing-masing memiliki orde 𝑛 × 𝑛. Bukti: 1. Diketahui tr 𝐴 + tr 𝐵 Akan dibuktikan bahwa tr(𝐴 + 𝐵) Bukti: 𝑛
𝑛
tr 𝐴 + tr 𝐵 =
𝑎𝑖𝑖 + 𝑖=1 𝑛
=
𝑏𝑖𝑖 𝑖=1
(𝑎𝑖𝑖 + 𝑏𝑖𝑖 ) 𝑖=1
= tr(𝐴 + 𝐵). 2. Diketahui 𝜆 tr(𝐴) Akan dibuktikan bahwa tr(𝜆𝐴) Bukti: 𝑛
𝜆 tr 𝐴 = 𝜆
𝑎𝑖𝑖 𝑖=1 𝑛
=
𝜆 𝑎𝑖𝑖 𝑖=1
= tr(𝜆𝐴). 3. Diketahui tr(𝐴) Akan dibuktikan bahwa tr(𝐴′) Bukti: 𝑛
tr 𝐴 =
𝑎𝑖𝑖 𝑖=1
= 𝑎11 + 𝑎22 + 𝑎33 + ⋯ + 𝑎𝑛𝑛 = tr(𝐴′). 4. Diketahui tr(𝐴𝐵) Akan dibuktikan bahwa tr(𝐵𝐴) Bukti: 𝑛
tr 𝐴𝐵 =
(𝐴𝐵)𝑖𝑖 𝑖=1 𝑛
=
(𝐴)𝑖 + (𝐵)𝑖 𝑖=1 𝑛 𝑛
=
𝑎𝑖𝑗 𝑏𝑗𝑖 . 𝑖=1 𝑗 =1
17
Lampiran 5 (Lanjutan)
𝑛
𝑛
=
𝑏𝑗𝑖 𝑎𝑖𝑗 𝑗 =1 𝑖=1 𝑛
=
𝐵𝑗 𝐴𝑗 𝑗 =1 𝑛
=
(𝐵𝐴)𝑗𝑗 𝑗 =1
= tr(𝐵𝐴). Bukti lengkap. ∎
18
Lampiran 6 Pembuktian Pendugaan Parameter 𝜷 dengan Algoritma EM (ExpectationMaximization) E-Step (Tahapan Expectation) Diketahui persamaan model regresi laten: 𝒀 = 𝑿𝜷 + 𝜺 Fungsi likelihood untuk persamaan model regresi: 𝐿 𝜷, 𝚿 = 2𝜋
−
𝑛𝑝 2
𝚿
𝑛 − 2 exp {−𝑡𝑟[𝚿−1
𝒀 − 𝑿𝜷 ′ (𝒀 − 𝑿𝜷)]/2}
Fungsi loglikelihood persamaan di atas menjadi: 𝑛𝑝
𝑛
ln𝐿 𝜷, 𝚿 = ln [ 2𝜋 − 2 𝚿 − 2 exp {−𝑡𝑟[𝚿−1 𝒀 − 𝑿𝜷 ′ (𝒀 − 𝑿𝜷)]/2}] 𝑛𝑝 𝑛 = − ln 2𝜋 − ln 𝚿 − 0.5𝑡𝑟[𝚿−1 𝒀 − 𝑿𝜷 ′ (𝒀 − 𝑿𝜷)] 2 2 𝚿−1 merupakan matriks kovarian simetri yang berorde 𝑝 × 𝑝, sehingga (𝚿 −1 )′ = (𝚿 −1 ). Ekspektasi dari fungsi loglikelihood bersyarat Y : 𝑛𝑝 𝑛 𝐸(ln𝐿(𝜷; 𝚿)|𝒀) = 𝐸{− ln 2𝜋 − ln 𝚿 − 0.5𝑡𝑟[𝚿−1 𝒀 − 𝑿𝜷 ′ (𝒀 − 𝑿𝜷)]|𝒀} 2 2 𝑛𝑝 𝑛 = − ln 2𝜋 − ln 𝚿 − 0.5𝐸{𝑡𝑟[𝚿−1 𝒀 − 𝑿𝜷 ′ (𝒀 − 𝑿𝜷)]|𝒀} 2 2 berdasarkan Lampiran 1, diperoleh: 𝑛 𝑛𝑝 ln 2𝜋 − ln 𝚿 − 0.5 𝑡𝑟 𝚿−1 {𝒀′ 𝒀 − 𝒀′ 𝑿𝜷 − 𝜷′𝑿′𝒀 + 𝜷′ (𝑿′𝑿)𝜷} 2 2 𝑛𝑝 𝑛 = − ln 2𝜋 − ln 𝚿 − 0.5 𝑡𝑟 𝚿−1 𝒀′ 𝒀 − 𝚿−1 𝒀′ 𝑿𝜷 − 𝚿 −1 𝜷′𝑿′𝒀 + 𝚿−1 𝜷′ (𝑿′𝑿)𝜷 2 2
=−
… (i)
berdasarkan Lampiran 2, diperoleh: 𝑛𝑝 𝑛 ln 2𝜋 − ln 𝚿 − 0.5 𝑡𝑟[𝚿−1 𝒀′ 𝒀 − 𝚿−1 𝜷′ 𝑿′ 𝑿 𝜷 − 𝚿−1 𝜷′𝑿′𝒀 + 𝚿−1 𝜷′ (𝑿′𝑿)𝜷] 2 2 𝑛𝑝 𝑛 = − ln 2𝜋 − ln 𝚿 − 0.5𝑡𝑟[𝚿−1 𝒀′ 𝒀 − 𝚿 −1 𝜷′ 𝑿′ 𝑿 𝜷 − 𝚿−1 𝜷′ 𝑿′ 𝑿 𝜷 + 𝚿−1 𝜷′ (𝑿′ 𝑿)𝜷] 2 2 …(ii) dari persamaan …(i)
=−
𝐸(ln𝐿(𝜷; 𝚿)|𝒀) 𝑛𝑝 𝑛 = − ln 2𝜋 − ln 𝚿 − 0.5 𝑡𝑟 𝚿−1 𝒀′ 𝒀 − 𝚿−1 𝒀′ 𝑿𝜷 − 𝚿 −1 𝜷′𝑿′𝒀 + 𝚿−1 𝜷′ (𝑿′𝑿)𝜷 2 2 berdasarkan Teorema trace, diperoleh: 𝑛𝑝 𝑛 ln 2𝜋 − ln 𝚿 − 0.5 𝑡𝑟[𝚿−1 𝒀′ 𝒀] + 0.5 𝑡𝑟[𝚿−1 𝒀′ 𝑿𝜷] + 0.5 𝑡𝑟[𝚿−1 𝜷′𝑿′𝒀] 2 2 −0.5 𝑡𝑟[𝚿−1 𝜷′ 𝑿′ 𝑿 𝜷] 𝑛𝑝 𝑛 = − ln 2𝜋 − ln 𝚿 − 0.5 𝑡𝑟[𝚿−1 𝒀′ 𝒀] + 0.5 𝑡𝑟[𝚿−1 𝒀′ 𝑿𝜷] + 0.5 𝑡𝑟[(𝚿−1 𝜷′𝑿′𝒀)′] 2 2 −0.5 𝑡𝑟[𝚿−1 𝜷′ 𝑿′ 𝑿 𝜷] 𝑛𝑝 𝑛 = − ln 2𝜋 − ln 𝚿 − 0.5 𝑡𝑟[𝚿−1 𝒀′ 𝒀] + 0.5 𝑡𝑟[𝚿−1 𝒀′ 𝑿𝜷] + 0.5 𝑡𝑟[𝒀′ 𝑿𝜷(𝚿−1 )′] 2 2 −0.5 𝑡𝑟[𝚿−1 𝜷′ 𝑿′ 𝑿 𝜷] =−
19
Lampiran 6 (Lanjutan) 𝑛𝑝 𝑛 ln 2𝜋 − ln 𝚿 − 0.5 𝑡𝑟[𝚿−1 𝒀′ 𝒀] + 0.5 𝑡𝑟[𝚿−1 𝒀′ 𝑿𝜷] + 0.5 𝑡𝑟[𝒀′ 𝑿𝜷𝚿−1 ] 2 2 −0.5 𝑡𝑟[𝚿−1 𝜷′ 𝑿′ 𝑿 𝜷] 𝑛𝑝 𝑛 = − ln 2𝜋 − ln 𝚿 − 0.5 𝑡𝑟[𝚿−1 𝒀′ 𝒀] + 0.5 𝑡𝑟[𝚿−1 𝒀′ 𝑿𝜷] + 0.5 𝑡𝑟[𝚿−1 𝒀′ 𝑿𝜷] 2 2 −0.5 𝑡𝑟[𝚿−1 𝜷′ 𝑿′ 𝑿 𝜷] 𝑛𝑝 𝑛 = − ln 2𝜋 − ln 𝚿 − 0.5 𝑡𝑟 𝚿−1 𝒀′ 𝒀 + 𝑡𝑟 𝚿 −1 𝒀′ 𝑿𝜷 − 0.5 𝑡𝑟[𝚿−1 𝜷′ 𝑿′ 𝑿 𝜷] 2 2
=−
M-Step (Tahapan Maximization) Pendugaan parameter 𝜷 diperoleh dengan memaksimumkan harapan dari fungsi loglikelihood bersyarat Y : 𝜕(𝐸(ln𝐿(𝜷; 𝚿)|𝒀)) =0 𝜕𝜷 𝑛𝑝 𝑛 𝜕(− 2 ln 2𝜋 − 2 ln 𝚿 − 0.5 𝑡𝑟 𝚿 −1 𝒀′ 𝒀 + 𝑡𝑟 𝚿 −1 𝒀′𝑿𝜷 − 0.5 𝑡𝑟[𝚿 −1 𝜷′ 𝑿′𝑿 𝜷]) ↔ =0 𝜕𝜷 𝜕(𝑡𝑟 𝚿 −1 𝒀′ 𝑿𝜷 − 0.5 𝑡𝑟[𝚿 −1 𝜷′ 𝑿′ 𝑿 𝜷] ↔0−0−0+ =0 𝜕𝜷 𝜕(0.5 𝑡𝑟[𝚿 −1 𝜷′ 𝑿′𝑿 𝜷]) ′ =0 ↔ 𝚿 −1 𝒀′ 𝑿 − 𝜕𝜷 𝜕(0.5 𝑡𝑟[𝚿 −1 𝜷′ 𝑿′ 𝑿 𝜷]) =0 ↔ 𝑿′𝒀(𝚿−1 )′ − 𝜕𝜷 −1 ′ 𝜕(0.5 𝑡𝑟[𝚿 𝜷 𝑿′𝑿 𝜷]) ↔ 𝑿′𝒀𝚿 −1 − =0 𝜕𝜷 Penjelasan: 𝑑[0.5 𝑡𝑟[𝚿 −1 𝜷′ 𝑿′ 𝑿 𝜷]] = 0.5 𝑡𝑟[𝚿 −1 (𝑑𝜷)′ 𝑿′ 𝑿 𝜷 + [𝚿 −1 𝜷′ 𝑿′𝑿 (𝒅𝜷)]] = 0.5 𝑡𝑟[𝚿 −1 (𝑑𝜷)′ 𝑿′ 𝑿 𝜷 + 0.5tr[𝚿 −1 𝜷′ 𝑿′ 𝑿 (𝒅𝜷)] = 0.5 𝑡𝑟[𝜷 𝑿′𝑿 (𝑑𝜷)′ (𝚿 −1 )′) + 0.5tr[𝚿 −1 𝜷′ 𝑿′ 𝑿 (𝒅𝜷)] = 0.5 𝑡𝑟[𝚿 −1 𝑿′ 𝑿 𝜷′(𝑑𝜷)] + 0.5 tr[𝚿−1 𝜷′ 𝑿′ 𝑿 (𝒅𝜷)] = 0.5 tr[𝚿 −1 𝜷′ 𝑿′ 𝑿 (𝒅𝜷)] + 0.5 tr[𝚿−1 𝜷′ 𝑿′𝑿 (𝒅𝜷)] = tr[𝚿 −1 𝜷′ 𝑿′𝑿 ]𝒅𝜷
𝜕[tr[𝚿−1 𝜷′ 𝑿′ 𝑿 𝜷]] 𝜕𝜷
= [𝚿−1 𝜷′ 𝑿′ 𝑿 ]′ ′
= 𝑿′ 𝑿 𝜷(𝚿−1 )′ = 𝑿′ 𝑿 𝜷(𝚿 −1 )
diperoleh ↔ 𝑿′𝒀𝚿 −1 − 𝑿′ 𝑿 𝜷(𝚿−1 ) = 0 ↔ 𝑿′ 𝒀𝚿−1 = 𝑿′ 𝑿 𝜷𝚿−1 ↔ 𝑿′ 𝑿 ↔ 𝑿′ 𝑿
−1 −1
𝑿′ 𝒀𝚿−1 𝚿 = 𝑿′ 𝑿 𝑿′ 𝒀 = 𝜷
−1
𝑿′ 𝑿 𝜷𝚿−1 𝚿
20
Lampiran 6 (Lanjutan) dari persamaan …(ii) 𝐸(ln𝐿(𝜷; 𝚿)|𝒀) 𝑛𝑝 𝑛 = − ln 2𝜋 − ln 𝚿 − 0.5𝑡𝑟[𝚿−1 𝒀′ 𝒀 − 𝚿 −1 𝜷′ 𝑿′ 𝑿 𝜷 − 𝚿 −1 𝜷′ 𝑿′ 𝑿 𝜷 + 𝚿−1 𝜷′ (𝑿′ 𝑿)𝜷] 2 2 berdasarkan Teorema trace, diperoleh: 𝑛𝑝 𝑛 ln 2𝜋 − ln 𝚿 − 0.5𝑡𝑟(𝚿 −1 𝒀′ 𝒀) + 0.5𝑡𝑟(𝚿−1 𝜷′ 𝑿′ 𝑿 𝜷) + 0.5𝑡𝑟(𝚿 −1 𝜷′ 𝑿′ 𝑿 𝜷) 2 2 −0.5𝑡𝑟(𝚿−1 𝜷′ (𝑿′ 𝑿)𝜷) 𝑛𝑝 𝑛 = − ln 2𝜋 − ln 𝚿 − 0.5𝑡𝑟(𝚿 −1 𝒀′ 𝒀) + 0.5𝑡𝑟(𝚿−1 𝜷′ 𝑿′ 𝑿 𝜷) + 0.5𝑡𝑟([𝚿−1 𝜷′ 𝑿′ 𝑿 𝜷]′) 2 2 −0.5𝑡𝑟(𝚿−1 𝜷′ (𝑿′ 𝑿)𝜷) 𝑛 𝑛𝑝 = − ln 2𝜋 − ln 𝚿 − 0.5𝑡𝑟(𝚿 −1 𝒀′ 𝒀) + 0.5𝑡𝑟(𝚿−1 𝜷′ 𝑿′ 𝑿 𝜷) + 0.5𝑡𝑟(𝚿 −1 𝜷′ 𝑿′ 𝑿 𝜷) 2 2 −0.5𝑡𝑟(𝚿−1 𝜷′ (𝑿′ 𝑿)𝜷) 𝑛𝑝 𝑛 = − ln 2𝜋 − ln 𝚿 − 0.5𝑡𝑟(𝚿 −1 𝒀′ 𝒀) + 𝑡𝑟(𝚿−1 𝜷′ 𝑿′ 𝑿 𝜷) − 0.5𝑡𝑟(𝚿−1 𝜷′ (𝑿′ 𝑿)𝜷) 2 2
=−
Pendugaan parameter 𝜷 diperoleh dengan memaksimumkan harapan dari fungsi loglikelihood bersyarat Y : 𝜕(𝐸(ln𝐿(𝜷; 𝚿)|𝒀)) =0 𝜕𝜷 ′ 𝑛𝑝 𝑛 𝜕(− 2 ln 2𝜋 ) 𝜕(2 ln 𝚿 ) 𝜕(0.5𝑡𝑟(𝚿−1 𝒀′ 𝒀)) 𝜕(𝑡𝑟(𝚿−1 𝜷 𝑿′ 𝑿 𝜷)) ↔ − − + 𝜕𝜷 𝜕𝜷 𝜕𝜷 𝜕𝜷 −1 𝜕(0.5𝑡𝑟(𝚿 𝜷′ (𝑿′ 𝑿)𝜷)) − =0 𝜕𝜷 ′ 𝜕(𝑡𝑟(𝚿−1 𝜷 𝑿′ 𝑿 𝜷)) 𝜕(0.5𝑡𝑟(𝚿−1 𝜷′ (𝑿′ 𝑿)𝜷)) ↔ − =0 𝜕𝜷 𝜕𝜷 Penjelasan: ′
′
𝑑(𝑡𝑟(𝚿−1 𝜷 𝑿′ 𝑿 𝜷)) = 𝑡𝑟(𝚿−1 𝜷 𝑿′ 𝑿 )𝑑𝜷 ′
𝜕(𝑡𝑟(𝚿−1 𝜷 𝑿′ 𝑿 𝜷)) 𝜕𝜷
′
= (𝚿−1 𝜷 𝑿′ 𝑿 )′ −1
= 𝑿′ 𝑿 ′𝜷(𝚿 )′ = 𝑿′ 𝑿 𝜷𝚿−1
𝑑[0.5 𝑡𝑟[𝚿 −1 𝜷′ 𝑿′ 𝑿 𝜷]] = 0.5 𝑡𝑟[𝚿 −1 (𝑑𝜷)′ 𝑿′ 𝑿 𝜷 + [𝚿 −1 𝜷′ 𝑿′𝑿 (𝒅𝜷)]] = 0.5 𝑡𝑟[𝚿 −1 (𝑑𝜷)′ 𝑿′𝑿 𝜷 + 0.5tr[𝚿−1 𝜷′ 𝑿′ 𝑿 (𝒅𝜷)] = 0.5 𝑡𝑟[𝜷 𝑿′𝑿 (𝑑𝜷)′ (𝚿 −1 )′) + 0.5tr[𝚿−1 𝜷′ 𝑿′ 𝑿 (𝒅𝜷)] = 0.5 𝑡𝑟[𝚿 −1 𝑿′ 𝑿 𝜷′(𝑑𝜷)] + 0.5 tr[𝚿−1 𝜷′ 𝑿′ 𝑿 (𝒅𝜷)] = 0.5 tr[𝚿 −1 𝜷′ 𝑿′𝑿 (𝒅𝜷)] + 0.5 tr[𝚿−1 𝜷′ 𝑿′ 𝑿 (𝒅𝜷)] = tr[𝚿−1 𝜷′ 𝑿′𝑿 ]𝒅𝜷
Lampiran 6 (Lanjutan)
21
𝜕[tr[𝚿−1 𝜷′ 𝑿′ 𝑿 𝜷]] 𝜕𝜷
= [𝚿−1 𝜷′ 𝑿′ 𝑿 ]′ ′
= 𝑿′ 𝑿 𝜷(𝚿−1 )′ = 𝑿′𝑿 𝜷(𝚿−1 )
diperoleh ↔ 𝑿′ 𝑿 𝜷𝚿−1 − 𝑿′ 𝑿 𝜷(𝚿−1 ) = 0 ↔ 𝑿′ 𝑿 𝜷𝚿−1 = 𝑿′ 𝑿 𝜷𝚿−1 ↔ 𝑿′ 𝑿
−1
𝑿′ 𝑿 𝜷𝚿−1 𝚿 = 𝑿′ 𝑿
↔ 𝜷 = 𝜷.
Jadi pendugaan parameter, yaitu: 𝜷 = 𝜷 = 𝑿′ 𝑿 Bukti lengkap. ∎
−1
𝑿′ 𝒀.
−1
𝑿′ 𝑿 𝜷𝚿 −1 𝚿
22
Lampiran 7 Kode R untuk Gambar 1 x <- seq(0.0, 1.0, 0.0055) Density <- dbeta(x, 0.5, 0.3) plot(x, Density, type="l")
23
Lampiran 8 Kode R sebagai Pendefinisian Awal (Tarpey & Petkova) > # procedure EMNORM > # procedure EMNORM > > # input: k Number of Mixture components ># y N-vector of observed data ># prior k-vector of initial prior probabilities ># mu k by p matrix of initial means ># var k array of initial p by p covariance matrices ># pool indicator for equality of variances: ># pool<-0: unequal variances; pool<-1: equal variances > > # output prior final initial probabilities ># mu final means ># var final variances ># loglik value of (complete data) log-likelihood ># nit number of iterations of EM-algorithm ># post posterior probability estimates > > emnorm <- function(k,y, prior, mu, var, pool) +{ + N <- dim(y)[1] # number of observations + p <- dim(y)[2] # dimension of data + eps <- .000000001 # convergence criterion + change <- 1 # initial test value for convergence + maxiter <- 2000 # maximum number of iterations + nit <- 0 # initialize iteration counter + #param <- c(prior, mu, var) # arrange all parameters in a vector + post <- matrix(0, N, k) # open matrix for posterior probabilities + while (change > eps && nit <= maxiter) # start iterations +{ + options(digits = 6) # format for the display of numerical results + muold <- mu # store old parameter values + varold <- var + priorold <- prior + + # Evaluate component densities and the log-likelihood + f <- matrix(0,N,k) + i <- 0 + for(i in 1:k-1) +{ + i <- i+1 + expo <- 0.5*diag(as.matrix(sweep(y,2,mu[i,]))%*%solve(var[,,i])%*%as.matrix(t(sweep(y,2,mu[i,])))) + w <- eigen(var[,,i]) + f[,i] <- prior[i,1]*(2*pi)^(-p*.5)/sqrt(exp(sum(log(w$values))))*exp(expo) +} + loglik <- sum(log(as.matrix(apply(f,1,sum)))) # evaluate log-likelihood function + + # Compute posterier probabilities for each component + i <- 0 + for(i in 1:k-1) +{ + i <- i+1 + post[,i] <- f[,i]/as.matrix(apply(f,1,sum)) +}
24
Lampiran 8 (Lanjutan) + post <- as.matrix(post) + + + # comment next line out if you don't want a protocol of the algorithm + #cat(format(c(nit,loglik), justify = "right"), fill = T) + + prior <- as.matrix(apply(post, 2, mean)) # M-step: prior probabilities + + mu <- as.matrix(t(post)) %*%as.matrix(y)/N # M-step: means + i <- 0 + for(i in 1:k-1) +{ + i <- i+1 + mu[i,] <- mu[i,]/prior[i,1] +} + + + varcommon <- matrix(0,p,p) + i <- 0 + for(i in 1:k-1) +{ + i <- i+1 + var[,,i] <- 1/(N*prior[i,])*t(y)%*%as.matrix(diag(post[,i]))%*%as.matrix(y)as.matrix(mu[i,])%*%mu[i,] # m-step: variances + varcommon <- varcommon+var[,,i]*prior[i,1] +} + if (pool ==1) +{ + var <- array(varcommon, c(p,p,k)) # M-step: common variance +} + change1 <- max(abs(prior - priorold)) # test value for convergence + change2 <- max(abs(mu-muold)) + change3 <- max(abs(var-varold)) + change <- max(change1,change2, change3) + + + nit <- nit + 1 # increase interation counter +} + if (nit >= maxiter) + cat("algorithm did not converge in ") # warning message if convergence + maxiter # is not reached + " iterations" + results <- list("prior" = prior, "mu" = mu, "var " = var, + "loglik" = loglik, "nit" = nit, "post" = post) + results +}
25
Lampiran 9 Kode R untuk Algoritma K-Means, Maximum Likelihood dan Algoritma EM (Tarpey & Petkova) > source("emnorm.r") > y <- read.table("improvement.dat",header=F) > y <- as.matrix(y) > hist(y,20, freq=FALSE, xlab="Difference", main="Prozac Difference Data",ylab="Relative Frequency") > n<- dim(y)[1] > plot(density(y)) > t <- seq(0,1, by = .005) > shapiro.test(y) > k <- 2 > prior <- matrix(1,k)/k > pool <- 0 > S <- cov(y) > p <- dim(y)[2] > kmean <- kmeans(y, k, iter.max = 100) > mu <- as.matrix(sort(kmean$centers)) > var <- array(S, c(p,p, k)) > results <- emnorm(k,y, prior, mu, var, pool) > muvar<- t(matrix(cbind(results$prior,results$mu,results$var),6,1)) > s <- sqrt(sum(kmean$withinss)/n) # initial value of sigma > bhat <- matrix(1,2,1) > bhat[1,1] <- results$mu[1] # Initial values for beta0 and beta1 estimates > bhat[2,1] <- -results$mu[1]+results$mu[2] > ab <- matrix(1,2,1) > nit <- 1 > maxiter <- 200 > while ( nit <= maxiter) +{ + nit <- nit+1 + xx <- matrix(0,n,4) + for(i in 1:n) +{ + yi = y[i,1] + fd<- function(x) {dnorm(yi, mean=(bhat[1,1] +bhat[2,1]*x), sd=s, log = FALSE)*x^(ab[1,1]1)*(1-x)^(ab[2,1]-1)} + fx<- function(x) {dnorm(yi, mean=(bhat[1,1]+bhat[2,1]*x), sd=s, log = FALSE)*x^ab[1,1]*(1x)^(ab[2,1]-1)} + fx2<- function(x) {dnorm(yi, mean=(bhat[1,1]+bhat[2,1]*x), sd=s, log = FALSE)*x^(ab[1,1]+1)*(1-x)^(ab[2,1]-1)} + flx<- function(x) {dnorm(yi, mean=(bhat[1,1]+bhat[2,1]*x), sd=s, log = FALSE)*log(x)*x^(ab[1,1]-1)*(1-x)^(ab[2,1]-1)} + flx2<- function(x){dnorm(yi, mean=(bhat[1,1]+bhat[2,1]*x), sd=s, log = FALSE)*log(1x)*x^(ab[1,1]-1)*(1-x)^(ab[2,1]-1)} + denom <- integrate(fd, 0, 1) + numx <- integrate(fx, 0,1) + numx2 <- integrate(fx2,0,1) + numlx <- integrate(flx, 0,1) + numlx2<- integrate(flx2,0,1) + xx[i,1] <- numx$value/denom$value + xx[i,2] <- numx2$value/denom$value + xx[i,3] <- numlx$value/denom$value + xx[i,4] <- numlx2$value/denom$value +} + xprimex <- matrix(1,2,2) + xprimex[1,1] <- n
26
Lampiran 9 (Lanjutan) + xprimex[1,2] <- sum(xx[,1]) + xprimex[2,1] <- xprimex[1,2] + xprimex[2,2] <- sum(xx[,2]) + x0 <- matrix(1,n,2) + x0[,2] <- xx[,1] + bhat <- solve(xprimex)%*%t(x0)%*%y + s<-sqrt((t(y)%*%y - t(bhat)%*%xprimex%*%bhat)/n) + xbar <- as.matrix(apply(xx,2,mean)) + m1 <- xbar[1,1] + m2 <- xbar[2,1] + ab[2,1] <- m1*(1-m1)^2/(m2-m1^2) + m1 -1 + ab[1,1] <- ab[2,1]*m1/(1-m1) + S <- matrix(0,2,1) + Info <- matrix(0,2,2) + eps <- 0.0000001 + change <- 100 + while (change > eps) +{ + abold <- ab + S[1,1] <- n*(digamma(ab[1,1]+ab[2,1])-digamma(ab[1,1]))+sum(xx[,3]) + S[2,1] <- n*(digamma(ab[1,1]+ab[2,1])-digamma(ab[2,1]))+sum(xx[,4]) + Info[1,1] <- n*(trigamma(ab[1,1]+ab[2,1])-trigamma(ab[1,1])) + Info[2,2] <- n*(trigamma(ab[1,1]+ab[2,1])-trigamma(ab[2,1])) + Info[2,1] <- n*(trigamma(ab[1,1]+ab[2,1])) + Info[1,2] <- Info[2,1] + Info<- -Info + ab <- ab + solve(Info)%*%S + change <- t(ab-abold)%*%(ab-abold) +} + kb <- gamma(ab[1,1]+ab[2,1])/(gamma(ab[1,1])*gamma(ab[2,1])) + loglike <- -n/2*log(2*pi)-.5*n*log(s^2) -.5/s^2*(t(y)%*%y 2*t(bhat)%*%t(x0)%*%y+t(bhat)%*%xprimex%*%bhat)+n*log(kb) +(ab[1,1]-1)*sum(xx[,3]) +(ab[2,1]-1)*sum(xx[,4]) + yhatlatent <- bhat[1,1] + bhat[2,1]*x0[,2] + yhatmix <- results$mu[1,1]*results$post[,1] + results$mu[2,1]*results$post[,2] + mseLatent <- t(y-yhatlatent)%*%(y-yhatlatent)/(n-5) + mseMix <- t(y - yhatmix)%*%(y-yhatmix)/(n-5) + abline(0,1) + cat("iteration ", nit, "\n") + cat("loglike =", loglike, "\n") + cat("Parameter Estimate", "\n") + cat("beta0: ", bhat[1,1], "\n") + cat("beta1: ", bhat[2,1], "\n") + cat("alpha0: ", ab[1,1], "\n") + cat("alpha1: ", ab[2,1], "\n") + cat("sigma: ", s, "\n") + cat("\n") + cat("\n") + bpfit<-gamma(ab[1,1]+ab[2,1])/(gamma(ab[1,1])*gamma(ab[2,1]))*t^(ab[1,1]-1)*(1t)^(ab[2,1]-1) + plot(t,bpfit, type='l', xlim=c(0,1),ylim=c(0,9), xlab="x") +} > yd=as.matrix(seq(-20,30, by=.001)) > fy2=NULL
27
Lampiran 9 (Lanjutan) > for (ij in 1:length(yd)) +{ + int2 <- function(x){dnorm(yd[ij,1], mean=(bhat[1,1]+bhat[2,1]*x),sd=s,log = FALSE)*gamma(ab[1,1]+ab[2,1])/(gamma(ab[1,1])*gamma(ab[2,1]))*x^(ab[1,1]-1)*(1x)^(ab[2,1]-1)} + fy12 <- integrate(int2,0,1)$value + fy2 = rbind(fy2, fy12) +} > emr <- t(as.matrix(muvar)) > fmix <- 1/sqrt(2*pi)*(emr[1,1]*exp(-.5*(ydemr[3,1])^2/emr[5,1])/sqrt(emr[5,1])+emr[2,1]*exp(-.5*(yd-emr[4,1])^2/emr[6,1])/sqrt(emr[6,1])) > sd=as.numeric(sqrt(var(y))) > ybar=apply(y,2,mean) > fnor <- 1/sqrt(2*pi*sd^2)*exp(-.5*(yd-ybar)^2/sd^2) > s2 <- s^2 +bhat[2,1]^2*ab[1,1]*ab[2,1]/((ab[1,1]+ab[2,1]+1)*(ab[1,1]+ab[2,1])^2) > ybar2=bhat[1,1]+bhat[2,1]*(ab[1,1]/(ab[1,1]+ab[2,1])) > s2 <- as.numeric(s2) > fnor2 <- 1/sqrt(2*pi*s2)*exp(-.5*(yd-ybar2)^2/s2) > par(mfrow=c(2,1)) > dy=density(y,bw=1.94) > hist(y,8, freq=FALSE, xlim=c(-10,30), ylim=c(0,0.12),ylab="Density", xlab="Difference",main="Baseline - (Week 1) HAM-D") > lines(dy, type='l', lwd=2 ) > lines(yd, fy2, lty=2, lwd=2, col=2) > lines(yd, fmix, lty=5, lwd=2, col=3) > bpfit<-gamma(ab[1,1]+ab[2,1])/(gamma(ab[1,1])*gamma(ab[2,1]))*t^(ab[1,1]-1)*(1t)^(ab[2,1]-1) > plot(t,bpfit, type='l', xlim=c(0,1),ylim=c(0,7),xlab="x",ylab="", main="Latent Beta Density", lwd=3)
28
Lampiran 10 Kode R untuk Plot Kuantil-Kuantil Data Depresi y <- read.table("improvement.dat",header=F) qqnorm(y[ ,1], xlab = "HAM-D Differences", ylab = "Latent regression Quantiles") qqline(y[ ,1], col = 2)
