6
2. ANALISIS DATA LONGITUDINAL Data longitudinal merupakan salah satu bentuk data berkorelasi. Pada data longitudinal, peubah respon diukur pada beberapa titik waktu untuk setiap subyek. Dalam studi longitudinal dimungkinkan untuk mempelajari perubahan respon antar waktu beserta faktor-faktor yang mempengaruhi perubahan tersebut, baik pada level populasi maupun level individu. Data longitudinal dicirikan oleh fakta bahwa pengamatan berulang dalam subyek yang sama cenderung berkorelasi (Zeger et al. 1988), sehingga modelmodel untuk analisis data longitudinal harus mengenali hubungan antara pengamatan berkala dalam subyek yang sama (Laird & Ware 1982). Korelasi antar pengamatan berulang dapat dimodelkan secara eksplisit (melalui pola matriks kovarian), maupun secara implisit (melalui pengaruh acak). Untuk memodelkan keheterogenan antar subyek, ada dua pendekatan dalam analisis data longitudinal (Zeger et al. 1988). Pertama dengan memodelkan keheterogenan secara eksplisit, dikenal sebagai pendekatan spesifik subyek, misalnya melalui model campuran dimana pengaruh spesifik subyek diasumsikan mengikuti suatu sebaran parametrik tertentu. Untuk data longitudinal kontinu, model linear campuran dari Laird dan Ware (1982) merupakan model yang sering digunakan. Kedua, respon rataan populasi dapat dimodelkan sebagai fungsi dari kovariat tanpa secara eksplisit memperhitungkan keheterogenan dari subyek ke subyek. Pendekatan ini dikenal sebagai model rataan populasi. Dalam pendekatan ini, matriks kovarian dari peubah respon secara langsung dimodelkan melalui struktur kovarian bagi galat intra-subyek. Model spesifik subyek dikenal juga sebagai model bersyarat, sedangkan model rataan populasi sering disebut model marginal (Pinheiro 2006). Perbedaan mendasar dari kedua model di atas adalah model spesifik subyek memungkinkan inferensi terhadap subyek tertentu, sedangkan pada model rataan populasi tidak. 2.1. Model Linear Campuran Model linier campuran untuk peubah respon kontinu bagi subyek ke-i (i=1,2,...,n) adalah sebagai berikut (Laird & Ware 1982): Yi
X i β Z i b i ε i , i 1,, n
7
Dalam hal ini
( Yi1 ,, Yimi )
Yi
adalah vektor peubah respon dari subyek ke-i, n =
total banyaknya subyek, dan mi = banyaknya deret data longitudinal dari subyek ke-i. Adapun Xi dan Zi adalah matriks rancangan masing-masing berdimensi mi x p dan mi x q yang bersesuaian dengan vektor pengaruh tetap pengaruh acak bi(qx1), sedangkan
i
dan vektor
(px1)
merupakan vektor galat intra-subyek
berdimensi mi x 1. Pada model di atas, diasumsikan bahwa q vektor pengaruh acak bi menyebar normal ganda dengan nilaitengah 0 dan matriks kovarian D, yakni bi ~ N(0,D). Demikian pula
i
~ N(0,Ri), serta bi dan
i
saling bebas. Hal ini berimplikasi
sebaran marginal dari Yi adalah normal dengan nilaitengah E(Yi) = Xi
dan
matriks kovarian Vi = ZiDZi + Ri. Matriks D dan Ri keduanya merupakan matriks simetrik definit positif. 2.2. Metode Pendugaan Parameter Metode pendugaan yang umum digunakan untuk menduga parameter dalam model linier campuran adalah metode kemungkinan maksimum (maximum likelihood/ML) atau metode kemungkinan maksimum berkendala (restricted maximum likelihood/REML). Jika
adalah parameter pengaruh tetap, dan
adalah parameter kovarian,
fungsi kepekatan peluang normal ganda bagi Yi, f(yi| , ) adalah: f ( y i | β, θ )
(2 )
ni 2
Vi
1 2
exp[
1 2
X i β) Vi 1 (y i
(y i
X i β)]
Fungsi kemungkinannya dapat dituliskan sebagai: L(β, θ)
i
f (y i | β, θ)
i
(2 )
ni 2
Vi
1 2
exp[
1 2
(y i
X i β) Vi 1 (y i
X i β)]
Dengan demikian fungsi log-kemungkinannya adalah: l (β, θ)
ln L(β, θ)
n 2
ln( 2 )
1 2
ln Vi i
1 2
(y i
X i β) Vi 1 (y i
X i β)
(2.1)
dan
secara
i
sedangkan n = ni adalah banyaknya amatan dalam gugus data. Walaupun memungkinkan untuk menduga parameter
simultan dengan memaksimumkan fungsi di atas, banyak algoritma komputasi yang menyederhanakan optimasi dengan cara menyembunyikan (profiling out) parameter
dari fungsi log-kemungkinan.
8
Dengan asumsi
, dan sebagai akibatnya Vi diketahui, maka fungsi log-
kemungkinan, l( , ) menjadi fungsi dari q(β)
1 2
saja, yaitu:
X i β) Vi 1 (y i
(y i
X i β)
i
Dengan generalized least squares (GLS), penduga bagi
dapat diperoleh
secara analitik sebagai berikut: 1
βˆ
1
X i Vi 1 y i
X i Vi X i i
(2.2)
i
Penduga di atas memiliki sifat statistik yang diinginkan, yakni merupakan penduga tak bias linier terbaik (BLUE) bagi . dan pengaruh tetap
Pendugaan parameter kovarian
dengan asumsi
tidak diketahui diuraikan sebagai berikut. Pertama untuk menduga
dibentuk
fungsi profil log-kemungkinan lML( ), yaitu dengan menggantikan parameter dalam persamaan (2.1) dengan penduganya pada persamaan (2.2), yaitu: l ML (θ)
n 2
ln( 2 )
1 2
ri Vi 1ri
1 2
ln Vi i
i 1
sedangkan
ri
yi
1
Xi
X i Vi 1 y i
X i Vi X i i
i
Pada umumnya pemaksimuman lML( ) terhadap linier, dengan kendala terhadap
merupakan optimasi tak
sedemikian sehingga persyaratan definit positif
bagi matriks D dan Ri terpenuhi. Nilai dugaan bagi
dapat diperoleh dengan cara
iterasi sampai konvergen. Setelah penduga ML bagi
diperoleh melalui proses iterasi, nilai βˆ dapat
dihitung tanpa iterasi dengan menggunakan persamaan (2.3) dan (2.4) sebagai berikut:
ˆ V i
ˆZ Zi D i
ˆ R i
(2.3)
1
βˆ
ˆ X Xi V i i
ˆ 1y Xi V i i
1
i
i
(2.4)
9
Karena Vi digantikan oleh penduganya Vˆ i , maka βˆ pada persamaan (2.4) dikatakan sebagai penduga tak bias linier terbaik empirik (Empirical Best Linear Unbiased Estimator/EBLUE) bagi . Ragam bagi βˆ merupakan matriks ragam-peragam berdimensi pxp, yaitu: 1
var(βˆ )
ˆ 1X Xi V i i i
Karena tidak mempertimbangkan hilangnya derajat bebas sebagai akibat menduga , maka penduga ML bagi
merupakan penduga yang berbias. Untuk
mengeliminasi bias ini dikembangkan bentuk alternatif dari metode ML yakni pendugaan REML. Penduga REML bagi
diperoleh berdasarkan optimasi fungsi log-
kemungkinan REML sebagai berikut: l REML (θ)
(n p) 2
ln( 2 )
1 2
ln Vi i
ri Vi 1ri
1 2
ln X i Vi 1 X i
1 2
i
i
Deskripsi dan pembandingan berbagai metode pendugaan dalam model linier campuran dapat dijumpai misalnya dalam Searle et al. 1992. Nilai prediksi bagi pengaruh acak merupakan nilai harapan bersyarat dari pengaruh acak jika nilai peubah respon diketahui, yang dapat dinyatakan sebagai: bˆ i
E(b i | Yi
yi )
ˆZ V ˆ 1 D i i (y i
X i βˆ )
Nilai harapan beryarat di atas merupakan prediktor tak bias linear terbaik empiris (Empirical Best Linear Unbiased Predictor/EBLUP) bagi pengaruh acak bi, karena diperoleh berdasarkan nilai dugaan matriks kovarian Vˆ i . Adapun matriks kovarian bagi prediktor pengaruh acak bˆ i adalah: Var (bˆ i )
ˆ Z (V ˆ D i i
1
ˆ 1X ( V i i
ˆ 1X ) 1 X V ˆ 1 ˆ Xi V i i i i )Z i D i
2.3. Pengujian Hipotesis dan Pembandingan Model Tersarang Hipotesis dari dua model yang memiliki hubungan tersarang dapat dibuat menjadi suatu formula. Model reference (model penuh) merupakan model yang lebih umum yang mencakup kedua hipotesis (H0 dan Ha), sedangkan model yang hanya mencakup H0 disebut model nested (model tersarang). Model penuh mengandung semua parameter yang diuji sedangkan model tersarang hanya
10
mengandung sebagian dari parameter tersebut. Uji yang digunakan untuk membandingkan kedua model tersebut adalah Likelihood Ratio Tests (LRTs). LRTs merupakan suatu uji yang membandingkan nilai fungsi likelihood untuk kedua model dengan persamaan: -2log(
Ltersarang L penuh
) = -2 log (Ltersarang) – (-2log(Lpenuh))~
2 df
sedangkan: Ltersarang = nilai fungsi likelihood pada model tersarang Lpenuh
= nilai fungsi likelihood pada model penuh
df
= selisih banyaknya parameter antara model penuh dan model tersarang LRTs juga dapat digunakan untuk menguji hipotesis parameter acak dan
tetap di dalam model. Pengujian parameter tetap dalam model menggunakan pendugaan ML, sedangkan dalam pengujian parameter acak digunakan pendugaan REML. Statistik ujinya adalah selisih (-2 ML/REML log likelihood) antara model penuh dan model tersarang seperti dinyatakan dalam persamaan di atas. 2.4. Penerapan terhadap Data Kasus HIV/AIDS Data yang digunakan untuk analisis data longitudinal merupakan hasil suatu percobaan klinis untuk membandingkan kemanjuran dan keamanan dua jenis obat antiretroviral dalam menangani pasien-pasien yang gagal atau tidak toleran terhadap terapi zidovudine (AZT). Percobaan melibatkan n = 467 pasien terinfeksi HIV yang terdiagnosa sebagai penderita AIDS atau memiliki jumlah sel CD4+ ≤ 300 per ml3 darah. Pasien dibagi secara acak untuk menerima salah satu dari dua jenis obat, yaitu didanosine (ddI) atau zalzitabine (ddC). Banyaknya sel CD4+ dicatat pada saat terlibat dalam studi (t = 0), dan kunjungan pada bulan ke 2, 6, 12 dan 18, sehingga maks mi = 5. Data ini digunakan oleh Guo dan Carlin (2004) untuk pemodelan bersama data longitudinal dan data daya tahan hidup (waktu sampai
terjadinya
kematian)
dari
penderita
HIV.
Data
diambil
dari
http://www.biostat.umn.edu/~brad/software.html. Peubah penjelasnya adalah Drug (ddI = 1, ddC = 0), Gender (male = 1, female = -1), PrevOI (AIDS diagnosis at study entry = 1, no AIDS diagnosis = -1), dan Stratum (AZT failure = 1, AZT intolerance = -1).
11
Sebelum dimodelkan dengan model linier campuran, terlebih dahulu dilakukan eksplorasi terhadap data. Boxplot banyaknya sel CD4+ pada lima titik waktu pengamatan untuk kedua jenis obat disajikan pada Gambar 2.1. Dari Gambar 2.1 tampak bahwa sebaran banyaknya sel CD4+ sangat menjulur ke kanan dengan banyak pencilan, mengindikasikan perlunya dilakukan transformasi data sebelum analisis berikutnya. Boxplot of CD4 0
2
ddC
600
6
12
18
ddI
500
CD4
400 300 200 100 0 0
2
6
12
18 Obstime
Panel variable: Drug-Type
Gambar 2.1. Boxplot data asal Transformasi akar dipilih untuk mengurangi kemenjuluran pola sebaran sekaligus untuk menstabilkan ragam, juga karena datanya merupakan data cacahan. Boxplot setelah data ditransformasi dapat dilihat pada Gambar 2.2. Setelah ditransformasi data terlihat lebih homogen serta lebih simetrik. Boxplot of Sqrt(CD4) 0 ddC
25
2
6
12
ddI
Sqrt(CD4)
20 15 10 5 0 0
2
6
12
18 Obstime
Panel variable: Drug Type
Gambar 2.2. Boxplot data hasil transformasi akar
18
12
Efek pengobatan umumnya tidak sama antar waktu, yaitu memungkinkan adanya interaksi antara jenis obat dengan waktu pengamatan. Pemeriksaan interaksi antara jenis obat dengan waktu pengamatan secara grafis disajikan melalui plot interaksi data hasil transformasi pada Gambar 2.3. Dari Gambar 2.3 dapat dilihat adanya perbedaan pola jumlah sel CD4+ antar waktu untuk kedua jenis obat. Untuk kelompok ddI, terjadi kenaikan jumlah sel CD4+ begitu diberikan obat ddI sampai bulan ke-2, namun turun lagi pada bulan ke-6, naik lagi sedikit pada bulan ke-12, kemudian terus menurun sampai bulan ke-18. Adapun untuk kelompok ddC terjadi penurunan jumlah sel CD4+ sampai bulan ke-6, namun kemudian jumlah sel CD4+ naik terus sampai bulan ke-18. Berdasarkan hasil ini efek interaksi akan dimasukkan dalam pemodelan. Interaction Plot for Sqrt(CD4) Data Means
8,0
Drug-Type ddC ddI
Mean
7,5
7,0
6,5
6,0 0
2
6 Obstime
12
18
Gambar 2.3. Plot interaksi antara waktu pengamatan dengan jenis obat Data longitudinal hasil transformasi akar banyaknya sel CD4+ dalam submodel-1 selanjutnya dimodelkan sebagai model linier campuran dengan persamaan sebagai berikut: wij
01
11Time ij
51 Stratum i
i 1,2, ,467
sedangkan b
b0 i
21Time ij
Drug i
b1i Timeij
31Genderi
41
Pr evOI
ij ,
j 1,2, , mi
(b0i , b1i ) ~ N 2 ( η, Σ) dan
Pada persamaan di atas, β (
ij
~ N (0,
01 ,
11 ,
21 ,
2
) 31 ,
41 ,
51 )
merupakan parameter
efek tetap, sedangkan b (b0i , b1i ) merupakan parameter efek acak untuk pasien ke-
13
i. Dalam hal ini b0i merupakan intersep acak untuk subyek ke-i, dan b1i adalah laju perubahan peubah respon per satuan waktu untuk pasien ke-i. Adapun
ij
merupakan galat intra-subyek yang diasumsikan menyebar normal dengan ragam yang sama. Hasil pemodelan dengan menggunakan model linier campuran disajikan pada Tabel 2.1, sedangkan output SAS disajikan pada Lampiran 1. Tabel 2.1. Nilai dugaan parameter beserta hasil uji dan SK 95% Parameter
Nilai dugaan
Galat baku
t
Intercept (β01) Time (β11) Time x Drug (β21)
8.0129 -0.1668 0.02998
0.3511 0.02038 0.02891
22.82 -8.19 1.04
<.0001 <.0001 0.3003
7.3230 -0.2069 -0.02682
8.7027 -0.1268 0.08678
Gender (β31)
-0.1582
0.3249
-0.49
0.6265
-0.7965
0.4800
PrevOI (β41)
-2.3152
0.2382
-9.72
<.0001
-2.7831
-1.8474
Stratum (β51)
-0.1309
0.2352
-0.56
0.5780
-0.5929
0.3311
bo
15.9111
1.1702
13.60
<.0001
13.8453
18.4789
σbo,b1
-0.1300
0.06169
-2.11
0.0350
-0.2509
-0.00913
0.02854
0.005968
4.78
<.0001
0.01969
0.04509
3.0716
0.1713
17.93
<.0001
2.7617
3.4370
σ σ
2
2 b1
σ2
Nilai-p
SK 95%
Berdasarkan Tabel 2.1 dapat dilihat bahwa peubah bebas yang berpengaruh nyata pada banyaknya sel CD4+ penderita HIV adalah obstime dan prevOI dengan nilai-p kurang dari 0.0001. Peubah prevOI yang nyata menunjukkan bahwa penderita yang terdeteksi AIDS pada awal studi memiliki jumlah sel CD4+ lebih rendah dibandingkan yang tidak terdeteksi AIDS, dengan rata-rata perbedaan jumlah sel CD4+ antara pasien yang tidak terdiagnosis AIDS pada awal studi dan yang terdeteksi AIDS sebesar 2.3152. Adapun peubah gender dan stratum tidak nyata pengaruhnya terhadap jumlah sel CD4+ pada α = 5%. Untuk kelompok obat ddI, nilai dugaan koefisien regresinya untuk Time sebesar -0.1668 + 0.02998 = -0.13682, sedangkan untuk kelompok ddC sebesar 0.1668. Dengan kata lain rata-rata penurunan jumlah sel CD4+ sebesar kelompok ddI sebesar 0.13682 per bulan, sedangkan untuk kelompok ddC sebesar 0.1668 per bulan. Namun perbedaan ini tidak nyata seperti dapat dilihat dari nilai-p sebesar 0.3003.
14
Semua komponen ragam pada model ini nyata pada taraf nyata 5%. Dari Tabel 2.1 diperoleh ragam jumlah sel CD4+ antar waktu untuk setiap pasien berkisar antara 2.7617 dan 3.4370 pada taraf kepercayaan 95% dengan nilai dugaan titik sebesar 3.0716. Nilai dugaan bagi ragam intersep sebesar 15.9111, dengan selang kepercayaan 95% yaitu (13.8453, 18.4789), yang berarti ada keragaman jumlah sel CD4+ awal antar pasien sewaktu masuk dalam studi. Ragam slope juga nyata dengan nilai dugaan ragam sebesar 0.02854, artinya laju penurunan jumlah sel CD4+ per bulan bervariasi antar pasien dengan keragaman berkisar antara 0.01969 dan 0.04509. Terdapat korelasi negatif antara intersep dan slope, yang ditunjukkan oleh nilai peragam antara intersep dan slope sebesar -0.13, atau korelasinya sebesar -0.193. Hasil pengujian nyata pada α = 5%, yang berarti penurunan jumlah sel CD4+ antar pasien dipengaruhi oleh jumlah sel CD4+ yang dimiliki sebelumnya (sewaktu masuk dalam studi). Semakin besar jumlah sel CD4+ awal yang dimiliki, semakin rendah laju penurunan jumlah sel CD4+ per bulan. Diagram pencar antara intersep dan slope serta boxplot untuk kedua pengaruh acak tersebut dapat dilihat pada Gambar 2.1. Marginal Plot of slope vs intersep
0,6
slope
0,4 0,2 0,0 -0,2 -0,4 -10
-5
0
5
10
15
intersep
Gambar 2.4. Korelasi antara intersep dan slope
