3. MODEL LINEAR CAMPURAN

3. MODEL LINEAR CAMPURAN 3.1. Sejarah Komponen Ragam Sejarah perkembangan komponen ragam diawali pada model dengan satu faktor. Formulasi pertama dari model efek acak dikemukakan oleh Airy pada tahun 1861. Airy melakukan pengamatan dengan teleskop pada fenomena yang sama untuk a malam dengan ni pengamatan pada malam-i. Airy menduga ragam 2 , dan menghitung rataannya galat untuk malam-i yang dinotasikan dengan σ ε,i

untuk untuk menduga σ ε2 . Fakta ini juga menunjukkan bahwa ahli astronomi pada masa itu telah mengenal konsep derajat bebas. Pengguna kedua model efek acak adalah Chauvenet pada tahun 1863 yang menurunkan formula untuk ragam rataan umum var( y..) , meskipun tanpa menuliskan modelnya (Searle et al, 1992). Perkembangan komponen ragam antara tahun 1900 sampai 1939 dipelopori oleh Fisher. Kontribusi besar Fisher (1925) pada model komponen ragam adalah metode analisis ragam (selanjutnya disebut metode ANOVA) yaitu menyamakan jumlah kuadrat dari analisis ragam dengan nilai harapannya (sesuai model acak atau campuran), sehingga menghasilkan satu gugus persamaan linear dalam komponen ragam yang akan diduga. Gagasan ini berasal dari pemakaian metode ANOVA untuk menurunkan penduga korelasi intra-kelas pada rancangan acak lengkap. Meskipun Fisher memperluas metode ANOVA satu faktor untuk data tak berimbang ke dalam model dua faktor dengan interaksi dan model-model yang kompleks, akan tetapi dalam papernya tidak pernah menggunakan operator E atau kalimat expected value untuk pendugaan komponen ragam (Searle et al, 1992). Fisher (1935) menggunakan istilah component of variation yang dirangkai dalam istilah error component pada paper-nya, yang tanpa diragukan merupakan asal pemakaian dari istilah variance component. Istilah komponen ragam pertama kali digunakan oleh Daniels (1939). Tippett (1931) menjelaskan metode ANOVA untuk pendugaan komponen ragam pada data berimbang dan perluasannya pada rancangan dua faktor tanpa interaksi pada model acak. Tippett mengabaikan kemungkinan efek interaksi pada model linear, sedangkan Fisher (1925) tidak hanya menggunakan istilah interaksi tetapi juga menjelaskan interaksi antara dua faktor A dan B. Tippett pertama kali

34 menampilkan tabel ANOVA dengan nilai harapan dari kuadrat tengah (Sahai dan Ageel, 2000). Komponen ragam untuk data tak berimbang pada model acak satu faktor dikemukakan oleh Cohran (1939). Perkembangan dari beragam model sangat intensif muncul antara 1940 sampai 1949. Eisenhart (1947) pertama kali menggunakan istilah Model I dan Model II untuk menyatakan model efek tetap dan model efek acak, dan juga memperhatikan kemungkinan dari model campuran. Sebelumnya Daniel (1939) dan Crump (1946) membahas tentang model acak, sedangkan Jackson (1939) membahas model campuran dua faktor tanpa interaksi. Crump juga menggunakan asumsi normal untuk pendugaan dengan metode ML. Kemudian diikuti Satterthwaite (1946) yang mengemukakan pendekatan pendugaan komponen ragam dengan sebaran contoh. Sedangkan Wald (1941) mengemukakan selang kepercayaan untuk ratio komponen ragam pada rancangan satu dan dua faktor pada data tak berimbang. Tahun 1950 sampai 1969 merupakan perkembangan terbesar dari metode pendugaan komponen ragam yang diawali dengan perluasan metode yg sudah ada dan diakhiri dengan pengukuhan metode baru berdasarkan ML atau

kriteria

norma minimum. Diawali Anderson dan Bancroft (1952) yang membahas komponen ragam pada data berimbang untuk model campuran dan model acak, data tak berimbang pada rancangan tersarang dan kelompok taklengkap. Dua tahun kemudian Bennett dan Franklin (1954) membahas nilai harapan kuadrat tengah dalam istilah komponen ragam. Henderson (1953) mengemukakan tentang pendugaan komponen ragam pada data tak berimbang. Paper ini dimotivasi oleh pekerjaan penulisnya yang berhubungan dengan analisis statistika pada perusahaan susu sapi. Paper ini mempresentasikan tiga metode yang berbeda dari penggunaan data tak berimbang pada model acak atau model campuran dengan sejumlah klasifikasi silang dan atau tersarang. Ketiga metode ini diadaptasi dari metode ANOVA dengan menyamaan (untuk data berimbang) jumlah kuadrat dengan nilai harapannya dan dikenal dengan Metode Henderson I, II, dan III. Pada tahun 1950-1969 banyak publikasi tentang sifat-sifat penduga dari metode ANOVA. Kelemahan penduga komponen ragam dengan metode ANOVA

35 pelan-pelan mulai disingkap. Hal ini tidak hanya karena metode Henderson I menghasilkan satu prosedur pendugaan komponen ragam untuk data tak berimbang yang sebelumnya tidak ada, tetapi juga ditunjang oleh kemajuan komputasi pada era awal komputer Diawali oleh Crump (1951) yang membahas metode ML untuk klasifikasi satu arah pada data berimbang dan tak berimbang, dan mendapatkan persamaan yang diselesaikan secara iteratif. Herbach (1959) menurunkan secara eksplisit penduga ML untuk model data berimbang. Kemudian diikuti Corbeil dan Searle (1976) yang menyarikan penduga ML pada sejumlah kasus data berimbang. Hartley dan Rao (1967) mengemukakan metode ML untuk model campuran dan model acak pada data berimbang atau tak berimbang. Kedatangan komputer dan perkembangan paket-paket komputasi sangat membantu dalam mengatasi permasalahan komputasi seperti prosedur iteratif metode ML. Miller (1973) meneliti sifat-sifat asimtotik dari penduga ML. Dengan demikian selang waktu antara penemuan metode ML oleh Fisher tahun 1925 dengan penerapannya pada pendugaan komponen ragam sekitar 40 tahun. Thompson (1962) juga mempelajari penduga ML dan mengenalkan gagasan maksimisasi pada bagian ML yang invarian terhadap parameter lokasi yaitu pada efek tetap dan dikenal dengan metode REML. Metode ini kemudian diterapkan untuk data tak berimbang oleh Patterson dan Thompson (1971). Kesulitan komputasi pada metode REML sama seperti pada metode ML, karena menurut Harville (1977) metode REML tidak lebih efektif dari metode ML. Hal yang menarik adalah solusi REML untuk data berimbang identik dengan penduga metode ANOVA. Berbeda dengan metode ML, metode REML memperhitungkan derajat bebas yang berhubungan dengan efek tetap. Metode lain dikemukakan oleh Rao (1972) yang dikenal dengan MINQUE (minimum norm quadratic unbiased). Metode ini tidak memerlukan asumsi sebaran bagi efek acak dan tidak memerlukan iterasi. Dari beberapa metode pendugaan komponen ragam, metode ML dan REML lebih disukai khususnya pada data berimbang.

36 3.2. Model Linear Campuran Bentuk umum model linear campuran adalah u    0  G y = X β + Z u + ε dengan   ~    ,    ε   0 0

0   R   

(1)

dimana X adalah matriks desain dari efek tetap yang teramati, β adalah vektor parameter pengaruh efek tetap yang tidak diketahui, Z adalah matriks desain efek acak yang teramati, u adalah vektor efek acak yang tidak diketahui, dan ε adalah vektor galat acak yang tidak diketahui. Sehingga nilaitengah dan matriks ragamperagam untuk y adalah E ( y ) = Xβ dan var (y ) = V = ZGZ' + R . Pendugaan Efek Tetap dan Efek Acak Persamaan model linear campuran didefinisikan sebagai :  X' R −1 X  −1  Z' R X

 β 

X' R −1 Z Z' R

−1

 X' R −1 y   −1   Z' R y 

   = Z + G −1  u 

(2)

Jika βˆ dan uˆ adalah solusi dari persamaan (2) maka baris kedua dari persamaaan

[

]

(2) adalah

Z' R − 1 Xβˆ + Z' R − 1 Z + G −1 uˆ = Z' R − 1 y

atau

uˆ = Z' R − 1 Z + G − 1

[

]

−1

(

Z' R −1 y − X βˆ

)

(3)

Sedangkan baris pertama dari persamaan (2) adalah X' R − 1 X βˆ + X' R − 1 Z uˆ = X' R −1 y

(4)

substitusi uˆ pada persamaan (4) menghasilkan

[

X' R −1 X βˆ + X' R −1 Z Z' R − 1 Z' + G −1

[

]

−1

X' R −1 X βˆ − X' R − 1 Z Z' R −1 Z' + G −1

atau

[

(

]

−1

= X' R −1 y − X' R − 1 Z Z' R −1 Z' + G − 1

atau

X' V -1Xβˆ = X' V −1y

sehingga

−1 βˆ = X' V -1X X' V −1y

(

)

Z' R −1 y − X βˆ = X' R −1 y

]

Z' R − 1 X βˆ

−1

Z' R −1 y

)

(5)

Penduga efek tetap βˆ adalah penduga GLS (generalized least squares) dan Xβˆ adalah BLUE (best linear unbiased estimator ) untuk Xβ . Bila y mempunyai

sebaran normal, maka βˆ adalah MLE (maximum likelihood estimator). Penduga efek acak uˆ pada persamaan (3) dapat dinyatakan sebagai

(

uˆ = GZ' V − 1 y − X βˆ

)

37 yang merupakan BLUP (best linear unbiased prediction) dari u (Christensen, 1987). Salah satu cara yang paling sederhana untuk mendapatkan BLUP adalah menggunakan

justifikasi

Henderson

(Henderson’s

justification)

dengan

menggunakan asumsi sebaran, yaitu y u ~ N ( Xβ + Zu, R ) dan

u ~ N(0, G)

maksimisasi fungsi kemungkinan bersama (y,u) pada β dan u yang tidak diketahui akan menghasilkan kriteria :

(y − Xβ − Zu )' R −1 (y − Xβ − Zu ) + u' G −1 u

(6)

Suku pertama pada persamaan (6) adalah jumlah kuadrat tertimbang dan suku keduanya adalah penalti. Hal ini menunjukkan bahwa kriteria BLUP dari (β,u) terdiri dari GLS ditambah dengan satu suku penalti. BLUP dari (β,u) dapat ditulis −1  βˆ  −1 −1   = C' R C + B C' R y ˆ u  

(

sebagai

)

dimana

0 C ≡ X Z dan B ≡  0

sehingga

BLUP (y ) = Xβˆ + Z uˆ = C C' R −1 C + B

[

0  G −1 

]

(

)−1 C' R −1y

Pendugaan Komponen Ragam Metode yang paling banyak digunakan dalam pendugaan matriks ragamperagam pada model linear campuran adalah metode ML dan REML, di samping metode lain seperti metode MINQUE dan

MIVQUE (minimum variance

quadratic unbiased). Metode REML menghasilkan penduga takbias bagi parameter matriks ragam-peragam, sedangkan ML menghasilkan penduga yang bias. Penduga ML dari V ditentukan berdasarkan model sebaran y ~ N ( Xβ + Zu , V )

Fungsi log-kemungkinan (log-likelihood) dari y berdasarkan model sebaran di atas adalah :

{

}

l ML (β, V ) = − 12 log V + (y − Xβ )' V −1 (y − Xβ ) + nlog (2 π )

MLE bagi (β , V )' diperoleh dengan

memaksimumkan l ML (β , V ) .

(7) Jika

dilakukan optimasi terhadap β dengan menganggap V tetap, maka akan diperoleh bahwa l ML (β, V ) maksimum pada

(

βˆ = X' V -1 X

)

−1

X' V −1y , yang sama dengan

38 BLUE pada persamaan (5). Substitusi βˆ pada persamaan (7) menghasilkan profil log-kemungkinan (profile log-likelihood) untuk V yaitu : ' l P ( V ) = − 12 log V + y − Xβˆ V −1 y − Xβˆ + n log (2 π )  

(

)

(

)

−1   = − 12 logV + y'V −1 Ι − X X'V −1X X'V−1 y − n2 log(2 π)    

[

]

(8)

Penduga ML bagi parameter-parameter dalam V dapat diperoleh dengan maksimisasi persamaan (8) terhadap parameter-parameternya. Prosedur iteratif pendugaan komponen ragam terdapat pada Searle et al (1992). Penurunan kriteria REML lebih rumit, yaitu meliputi maksimisasi fungsi kemungkinan dari kombinasi linear elemen y yang tidak tergantung pada β. Uraian lengkap tentang metode REML dapat dijumpai pada Christensen (1987) dan Searle et al (1992). Fungsi kriteria dari REML adalah ' l REML(V) = − 12 logV + log X'V−1X + y − Xβˆ V−1 y − Xβˆ + (n− p)log(2 π)  

(

) (

)

p

= l P (V ) − 12 log X' V −1X + 2 log(2 π )

Keuntungan utama dari REML dibanding ML adalah REML memperhitungkan derajat bebas dari efek tetap dalam model linear campuran. 3.3 Pemulusan dengan Model Linear Campuran Pada awalnya model linear campuran khususnya kasus model komponen ragam sukses diterapkan pada percobaan pemuliaan hewan yaitu untuk mengidentifikasi hewan yang mempunyai genetika unggul. Pada saat ini penerapan dalam bidang epidemiologi genetika lebih dikhususkan untuk tujuan tertentu pemuliaan hewan. Model linear campuran banyak diterapkan untuk analisa data longitudinal atau pengukuran berulang dimana galat acak berkorelasi (Milliken dan Johnson, 1984). Suatu paradigma baru di bidang pemulusan adalah pendekatan model linear campuran untuk pemulusan (Wand, 2003). Penerapan model linear campuran untuk metode pemulusan didorong oleh munculnya teori P-spline yang dikemukakan oleh Eiler dan Mark (1996) dan Ruppert dan Carroll (1997). Di samping itu, pendekatan ini juga ditunjang oleh munculnya paket program model linear campuran di SAS dan S-PLUS (Pedan, 2003).

39 P-spline adalah regresi spline yang diduga dengan minimisasi jumlah kuadrat terpenalti. Misalkan hubungan fungsional antara peubah penjelas x dengan respon y dimodelkan sebagai y i = s (xi ) + ε i

dan fungsi mulus s dimodelkan sebagai regresi spline yang dinotasikan sebagai s = X β* + Z u

dimana ε i bebas stokastik dengan nilaitengah nol dan ragam σ2, sedangkan

(

β * = β 0 , β1 ,L , β p

1  X = M 1 

(w )+p

x1

L

M

O

x1

L

)' dan u = (u p 1 ,..., u pK )'

 ( x1 − κ1)p+ x1p   dan  Z= M  M p  (x − κ )p xn  n 1+ 

adalah vektor parameter,

L ( x1 − κK )p+   adalah matriks desain, O M  L ( xn − κK )p+  

= w p Ι(w ≥ 0 ) adalah basis fungsi pangkat terpotong berderajat-p (truncated

power function, dan selanjutnya disingkat dengan FPT) dengan p ≥ 1 adalah bilangan bulat positif, dan κ 1 < ... < κ K adalah simpul tetap. Jumlah kuadrat terpenalti pada P-spline didefinisikan sebagai

( y − X β − Z u ) ( y − X β − Z u )+ λ u u '

*

'

*

(9)

Suku pertama pada persamaan (9) adalah jumlah kuadrat galat dan suku keduanya adalah penalti kekasaran, dan λ adalah parameter pemulus. Bila kriteria spline 2

terpenalti dibagi dengan σ ε2 dan σ u2 = σ ε , maka akan diperoleh λ 1

σ ε2

( y − X β − Z u ) ( y − X β − Z u )+ σ1 *

'

*

2 u

u 'u

(10)

Kriteria P-spline pada persamaan (10) sama dengan kriteria BLUP dari model linear campuran pada persamaan (6) bila matriks R = σ u2 I dan G = σ ε2 I . Efek samping yang menarik dari hubungan antara P-spline dengan model linear campuran adalah parameter pemulus P-spline merupakan rasio antara dua komponen ragam. Hubungan matematis antara P-spline dengan model linear campuran telah dibahas oleh beberapa peneliti antara lain Wang (1998), Fan dan Zhang (1998), Brumback et al (1999), Vebyla et al (1999), French et al (2001), dan Wand (2003).

40 Model linear campuran juga digunakan untuk prediksi spatial seperti yang dikemukakan oleh Diggle et al (1998) yaitu dengan mendesain komponen efek acak mempunyai fungsi ragam spatial (variogram) tertentu. Akhir-akhir ini juga banyak penelitian tentang prediksi spatial dimensi rendah untuk kriging dan spline-2 (thin-plate spline) dengan menggunakan model campuran. Pendekatan ini merupakan perluasan dari P-spline untuk data spatial yang dikemukakan oleh French et al (2001), Kamman dan Wand (2003), Ranalli et al (2005). Penerapan lainnya dari model linear campuran adalah untuk pendugaan area kecil (small area estimation) seperti yang dikemukakan oleh Kleinschmidt et al (2001) dan Opsomer (2005). Pada metode deret waktu, Tsimikas dan Ledolter (1997) merepresentasikan state-space dalam bentuk model linear campuran sehingga bisa diduga dengan metode REML. Pada saat ini banyak penelitian tentang pengembangan model linear campuran dengan menggabungkan pendekatan-pendekatan di atas seperti yang dikemukaan oleh Zhang et al (1997), Guo (2002), Wu dan Zhang (2002). Daftar Pustaka Anderson RL, Bancroft TL.1952. Statistical Theory in Research. New York : McGraw-Hill. Bennett CA, Franklin NL. 1954. Statistical Analysis in Chemistry and Chemical Industry. New York : John Wiley & Sons. Brumback BA, Ruppert D, Wand MP. 1999. Comment on Variable selection and function estimation in additive nonparametric regression using a data-based prior by Shively, Kohn and Wood. J Amer Stat Ass 94:794-797. Christensen R. 1984. Plane Answers to Complex Questions. The Theory of Linear Models. New York : Springer-Verlag. Cohran WG. 1939. The use of analysis of variance in enumeration by sampling. J Amer Stat Ass 34: 492-510. Corbeil RR, Searle SR. 1976. A comparison of variance component estimators. Biometrics 32: 779-791. Crump SL. 1946. The estimation of variance components analysis of variance. Biometrics Bull 2: 7-11. Crump SL. 1951. The variance component analysis. Biometrics 7 : 1-16. Daniel HE. 1939. The estimation of component of variance. J R Stat Soc Suppl 6:186-197.

41 Diggle PJ, Tawn JA, Moyeed RA. 1998. Model-based geostatistics (with discussion). Appl Stat 47:299-350. Eilers PHC, Marx BD. 1996. Flexible smoothing with B-splines and penalties (with discussion). Stat Sci 11: 89-121. Eisenhart C. 1947. The assumptions underlying the analysis of variance. Biometrics 3:1-21. Fan J, Zhang JT. 1998. Comment on Smoothing spline models for the analysis of nested and crossed samples of curves by Brumback and Rice. J Amer Stat Ass 93: 961-994. Fisher RA.1925. Statistical Methods for Research Workers. 1st edn. Edinburg and London: Olyver & Boyd. Fisher RA.1935. Discussion of Neyman et al. J R Stat Soc, Series B 2:154-155. French JL, Kammann EE, Wand MP. 2001. Comment on Semiparametric nonlinear mixed-effects models and their applications by Ke and Wang. J Amer Stat Ass 96:1285-1288. Guo W. 2002. Functional mixed effects models. Biometrics, 58, 121-128. Hartley HO, Rao JNK. 1967. Maximum likelihood estimation for the mixed analysis of variance model. Biometrika 54:93-108. Harville DA.1977. Maximum likelihood approaches to variance component estimation and to related problems. J Amer Stat Ass 72:320-340. Henderson CR. 1953. Estimation of variance and covariance components. Biometrics 9:226-252. Herbach LH. 1959. Properties of Model II type of analysis of variance tests. An optimum nature of F-test for Model II in balanced case. Ann Math Stat 30:030-959. Jackson RWB. 1939 Reliability of mental tests. Brit. J Psycol 29:267-287. Kammann EE, Wand MP. 2003. Geoadditive models. Appl Stat 52:1-18. Kleinschmidt I, Sharp BL, Clarke GPY, Curtis B, Fraser C. 2001. Use of Generalized linear mixed models in the spatial analysis of small-area malaria incidence rates in KwaZulu Natal, South Africa. Amer J Epidemiology 153(12): 1213-1221 Miller JJ. 1973. Asymptotic properties of maximum likelihood estimates in the mixed model of the analysis of variance. Ann Stat 5:746-762. Milliken GA, Johnson DE. 1984. Analysis of Messy Data, Volume I : Design Experiment. New York : Van Nostrand Reinhold Company. Opsomer JD. 2005. Small area estimation using penalized spline regression. International Biometric Society, Eastern North American Region meeting, Austin, TX, March 21, 2005. [terhubung berkala]. http://www.stat. colostate.edu/~nsu/starmap/pps/Opsomer.Spline_survey ENAR.pdf Pedan A. 2003. Smoothing with SAS PROC MIXED. SUGI 28: 268-28

42 Ranalli, MG, Breidt FJ, Wang H. 2005. Low-rank smoothing splines on complex domains. Seminar, Atlantic Ecology Division, EPA, Narragansett, RI. March 1, 2005. [terhubung berkala]. http://www.stat.colostate.edu/~nsu/ starmap/ pps/ ranalli_epaaed_march1.pdf Rao CR. 1972. Estimation of variance and covariance component in linear models. J Amer Stat Ass 67:112-115. Ruppert D, Carroll RJ. 1997. Penalized regression splines. Unpublished manuscript. [terhubung berkala]. http://www.orie.cornell.edu/~davidr/ papers/Index /index.html Sahai H, Ageel, MI. 2000. The Analysis of Variance. Fixed, Random, and Mixed Models. Boston : Birkhäuser. Satterthwaite FE. 1946. An approximate distribution of estimates of variance components. Biometric Bull 2:110-114. Searle SR, Casella G, McCulloch CE. 1992. Variance Component. New York : John Wiley & Sons. Thompson WA. 1962. The problem of components. Ann Math Stat 33:273-289.

negative estimates of variance

Tippett LHC. 1931. The Method of Statistics. 1st edn. London : Williams and Norgate. Tsimikas JF, Ledolter J. 1997. Mixed model representation of state space models : New smoothing results and their application to REML estimation. Stat Sinica 7: 973-991 Verbyla AP, Cullis BR, Kenward, MG, Welham SJ. 1999. The analysis of designed experiments and longitudinal data by using smoothing splines (with discussion). J R Stat Soc, Series C 48: 269-312. Wald A. 1941. On the analysis of variance in case of multiple classifications with unequal class frequencies. Ann Mathc Stat 12:346-350. Wand M. 2003. Smoothing and mixed models. Comp Stat 18 : 223–249. Wang Y, 1998. Mixed effects smoothing spline analysis of variance. J R Stat Soc, Series B 60:159-174. Wu H, Zhang JT. 2002 Local polynomial mixed-effects models for longitudinal data. J Amer Stat Ass 97:883-897. Zhang, D, Lin X, Raz, J, Sower MF. 1998. Semiparametric stochastic mixed models for longitudinal data. J Amer Stat Ass 93:710-719.