FItting Linear Mixed Model dengan metan
Shahrizal Muhammad Abdillah
Model Linier Campuran
Model linier campuran atau linier mixed model (LMM) adalah model perluasan dari regresi linier yang melibatkan efek tetap dan efek acak pada faktor-faktornya secara bersamaan di dalam model. Hal ini dilakukan karena secara alamiah, variabilitas yang muncul seringkali tidak bisa dijelaskan hanya dengan efek tetap. Selain itu juga, adanya struktur data yang hierarki, seperti perlakuan yang tersarang atau nested. Dengan kata lain, data memiliki tingkatan atau kelompok-kelompok yang saling terkait, sehingga tidak bisa dijelaskan hanya dengan efek tetap. Contohnya pada percobaan split plot, dimana pada petak utama (main plot) sering kali dianggap sebagai faktor acak, karena pada petak utama seringkali variasi yang ada muncul.
Percobaan MET
Percobaan MET atau multi-environment trials merupakan serangkaian percobaan yang melibatkan beberapa percobaan secara bersamaan di beberapa lingkungan yang diuji. Beberapa lingkungan tersebut memiliki perbedaan karakteristik kondisinya seperti iklim, jenis tanah, ketinggian, curah hujan, dll. Percobaan MET ini dilakukan untuk menduga interaksi GxE, yaitu varietas atau genotipe-genotipe uji ditanam di sejumlah lingkungan. Adanya interaksi GxE dapat mengakibatkan performa genotipe di setiap lingkungan berbeda-beda ataupun memiliki daya adaptasi yang rendah. Seiring dengan perubahan kondisi di lingkungan memungkinkan adaptabilitas genotipe menjadi rendah. Dengan mengasumsikan lingkungan sebagai efek acak pada percobaan, maka dalam model campuran ini, genotipe sebagai efek tetap.
BLUP
Pendugaan pada pengaruh acak disebut dengan prediksi, sedangkan pendugaan pada pengaruh tetap adalah estimasi. Dengan demikian prosedur yang berhubungan dengan prediksi disebut BLUP (best linier unbiased prediction) atau prediksi tak bias.
Dalam model campuran, respons (\(y\)) dinyatakan sebagai:
\(y=Xβ+Zu+ε\mathbf{y} = \mathbf{X} \boldsymbol{\beta} + \mathbf{Z} \mathbf{u} + \boldsymbol{\varepsilon}y=Xβ+Zu+ε\)
di mana:
\(Xβ\) = efek tetap (misalnya, rata-rata genotipe)
\(Zu\) = efek acak, dengan u∼N(0,G)
ε∼N(0,R) = error/residual
BLUP memprediksi u dengan menggunakan kombinasi data yang diamati (y) dan informasi dari struktur varians-kovarians efek acak (G dan R).
Formula BLUP untuk efek acak u:
\(\hat u = GZ^TV^{−1}(y−X \hatβ)\)
di mana:
\(V = ZGZ^T + R\) adalah matriks kovariansi total
G adalah varians efek acak
R adalah varians residual
Z adalah matriks desain untuk efek acak
Paket metan
Pada paket metan menyediakan analisis model linier
campuran baik pada percobaan MET ataupun percobaan tunggal.
Metode yang digunakan untuk menampilkan hasil analisis model linier
campuran pada percobaan MET adalah gamem_met().
untuk menampilkan plot residual dari masing-masing efek acak dapat
dengan plot(x, type="re") dimana x adalah objek dari hasil
eksekusi gamem_met() dengan re sebagai
residual. get_model_data() digunakan untuk menampilkan
beberapa hasil variabel yang telah dihitung dalam model linier
campuran.
Dataset
Pada dataset yang digunakan pada analisis ini adalah
data("dasilva.maize") yang tersedia pada paket
agridat. Dataset ini merupakan percobaan MET
jagung di texas dengan beberapa peubah respon pengamatan. Peubah respon
yang digunakan pada analisis ini adalah yield. Untuk lebih
jauh mengetahui dataset ini dapat memanggilnya dengan perintah
?dasilva.maize.
Analisis
Berikut adalah penerapan langsung coding R untuk model linier
campuran pada dataset dasilva.maize.
# Set up library dan data
library(agridat)
library(metan)## Registered S3 method overwritten by 'GGally':
## method from
## +.gg ggplot2## |=========================================================|## | Multi-Environment Trial Analysis (metan) v1.18.0 |## | Author: Tiago Olivoto |## | Type 'citation('metan')' to know how to cite metan |## | Type 'vignette('metan_start')' for a short tutorial |## | Visit 'https://bit.ly/pkgmetan' for a complete tutorial |## |=========================================================|data("dasilva.maize")
dat <- dasilva.maize
# Cek struktur data
str(dat)## 'data.frame': 1485 obs. of 4 variables:
## $ env : Factor w/ 9 levels "E1","E2","E3",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ rep : Factor w/ 27 levels "R1","R10","R11",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ gen : Factor w/ 55 levels "G01","G02","G03",..: 1 2 3 4 5 6 7 8 9 10 ...
## $ yield: num 7.91 6.41 5.86 7.77 3.96 ...# Fit model LMM
fit <- gamem_met(dat, env = env,
gen = gen,
rep = rep,
random="env",
resp=yield
)## Evaluating trait yield |=========================================| 100% 00:00:06 ## Method: REML/BLUP## Random effects: REP(ENV), ENV, GEN:ENV## Fixed effects: GEN## Denominador DF: Satterthwaite's method## ---------------------------------------------------------------------------
## P-values for Likelihood Ratio Test of the analyzed traits
## ---------------------------------------------------------------------------
## model yield
## COMPLETE NA
## REP(ENV) 3.69e-04
## ENV 4.67e-18
## GEN:ENV 3.87e-15
## ---------------------------------------------------------------------------
## All variables with significant (p < 0.05) genotype-vs-environment interaction# Tampilkan hasil model
fit## Variable yield
## ---------------------------------------------------------------------------
## Individual fixed-model analysis of variance
## ---------------------------------------------------------------------------
## NULL
## ---------------------------------------------------------------------------
## Fixed effects
## ---------------------------------------------------------------------------
## # A tibble: 1 × 7
## SOURCE `Sum Sq` `Mean Sq` NumDF DenDF `F value` `Pr(>F)`
## <chr> <dbl> <dbl> <int> <dbl> <dbl> <dbl>
## 1 GEN 318.7 5.902 54 432.0 5.061 1.478e-22
## ---------------------------------------------------------------------------
## Random effects
## ---------------------------------------------------------------------------
## # A tibble: 4 × 3
## Group Variance Percent
## <chr> <dbl> <dbl>
## 1 ENV 6.788 81.53
## 2 GEN:ENV 0.3351 4.025
## 3 REP:ENV 0.03687 0.4428
## 4 Residual 1.166 14.01
## ---------------------------------------------------------------------------
## Likelihood ratio test
## ---------------------------------------------------------------------------
## model npar logLik AIC LRT Df Pr(>Chisq)
## <none> 1 59 -2400.5 4918.9
## (1 | REP:ENV) 4 58 -2406.8 4929.6 12.681 1 0.0003694 ***
## (1 | ENV) 2 58 -2438.0 4991.9 75.014 1 < 2.2e-16 ***
## (1 | GEN:ENV) 3 58 -2431.3 4978.7 61.763 1 3.873e-15 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## ---------------------------------------------------------------------------
## Variance components and genetic parameters
## ---------------------------------------------------------------------------
## NULL
## ---------------------------------------------------------------------------
## Principal component analysis of the G x E interaction matrix
## ---------------------------------------------------------------------------
## NULL
## ---------------------------------------------------------------------------
## Some information regarding the analysis
## ---------------------------------------------------------------------------
## # A tibble: 12 × 2
## Parameters Values
## <chr> <chr>
## 1 Mean "7.75"
## 2 SE "0.07"
## 3 SD "2.81"
## 4 CV "36.25"
## 5 Min "2.35 (G19 in E5)"
## 6 Max "18.46 (G31 in E6)"
## 7 MinENV "E2 (4.55)"
## 8 MaxENV "E6 (13.19)"
## 9 MinGEN "G19 (6.1) "
## 10 MaxGEN "G06 (9.25) "
## 11 Ngen "55"
## 12 Nenv "9"Pada output fit diatas menunjukkan pendekatan untuk
menduga efek acak pada model ini adalah REML sedangkan metodenya adalah
BLUP. Efek acak yang terbaca adalah REP(ENV)
merepresentasikan ulangan dalam lingkungan (menunjukkan struktur
hierarkis data). ENV adalah lingkungan dan
GEN:ENV adalah interaksi Genotipe dengan Lingkungan.
Output hasil analisis model LMM menunjukkan seluruh komponen acak
yang untuk tes rasio Likelihood < 0.0001. Untuk F hitung
dari faktor tetap genotipe (GEN) sebesar 5.061, dengan DB
genotipe sebesar 54 sedangkan DB penyebut sebesar 432. Hasil ini
menunjukkan nilai P-value < 0.05 yang artinya ada perbedaan rataan
hasil yang signifikan antar genotipe yang diuji.
Selanjutnya pendugaan ragam pada komponen acak sebesar 6.788
ditunjukkan oleh lingkungan (ENV) dengan persentase sebesar
81.53% sedangkan interaksi GxE sebesar 0.3351 dengan persentase sebesar
4.025. Artinya variabilitas lingkungan memberikan pengaruh keragaman
tertinggi pada peubah respon hasil (yield) yang diukur
dibandingkan dengan interaksi GxE.
Hasil output LMM yang digunakan untuk mengevaluasi berbagai struktur efek acak menggunakan likelihood ratio test (LRT). Mari kita bahas setiap bagian secara sistematis.
Bagian-Bagian Output:
model: Menunjukkan model yang diuji, dengan berbagai efek acak yang dimasukkan.
npar (number of parameters): Jumlah parameter yang diestimasi dalam setiap model.
logLik (log-likelihood): Ukuran seberapa baik model cocok dengan data. Nilai lebih tinggi (kurang negatif) menunjukkan model yang lebih baik.
AIC (Akaike Information Criterion): Digunakan untuk membandingkan model; model dengan AIC lebih kecil lebih disukai.
LRT (Likelihood Ratio Test): Uji statistik yang membandingkan model dengan dan tanpa efek acak tertentu.
Df (Degrees of Freedom): Derajat kebebasan dalam uji LRT.
Pr(>Chisq): p-value untuk uji LRT. Nilai kecil (biasanya <0.05) menunjukkan bahwa efek acak yang diuji berkontribusi signifikan terhadap model.
Interpretasi Model
Model dasar ().Ini adalah model tanpa efek acak, hanya memiliki satu parameter yang diestimasi. Log-likelihood = -2400.5, AIC = 4918.9.
(1 | REP:ENV). Efek acak dari blok REP (Replicate) dalam ENV (Environment). Menambahkan efek acak REP dalam ENV menurunkan log-likelihood menjadi -2406.8, meningkatkan AIC menjadi 4929.6. LRT = 12.681, p-value = 0.0003694 (*, signifikan)** → Efek acak REP dalam ENV berpengaruh secara signifikan.
(1 | ENV) → Efek acak dari Environment (Lokasi/ENV). Menambahkan efek acak ENV menurunkan log-likelihood menjadi -2438.0, meningkatkan AIC menjadi 4991.9. LRT = 75.014, p-value = <2.2e-16 (*, sangat signifikan)** → Efek acak lingkungan memiliki pengaruh yang sangat signifikan terhadap model.
(1 | GEN:ENV) → Efek acak dari Interaksi Genotipe × Environment (GEN:ENV). Menambahkan efek acak GEN × ENV menurunkan log-likelihood menjadi -2431.3, meningkatkan AIC menjadi 4978.7. LRT = 61.763, p-value = 3.873e-15 (*, sangat signifikan)** → Interaksi Genotipe × Lingkungan berpengaruh signifikan terhadap model.
Kesimpulan
Semua efek acak yang diuji (REP:ENV, ENV, GEN:ENV) berkontribusi signifikan terhadap model (p-value < 0.05).
Efek lingkungan (ENV) memiliki dampak paling besar (LRT = 75.014, p-value < 2.2e-16), menunjukkan bahwa variasi antar lokasi sangat mempengaruhi hasil.
Efek interaksi GEN:ENV juga penting (LRT = 61.763, p-value = 3.873e-15), menunjukkan bahwa respons genotipe berbeda di berbagai lingkungan.
Efek blok dalam lingkungan (REP:ENV) juga signifikan, tetapi dampaknya lebih kecil dibanding ENV dan GEN:ENV.
Model terbaik bisa ditentukan dengan mempertimbangkan AIC yang lebih kecil dan log-likelihood yang lebih besar (kurang negatif).
Untuk menampilkan koefisien blup dari efek acak, yang salah satunya direpresentasikan oleh interaksi GxE dapat ditampilkan sebagai berikut.
BLUP <- waasb(dat,
resp = yield,
gen = gen,
env = env,
rep = rep)## Missing cells for: ENVE2:REPR1, ENVE3:REPR1, ENVE4:REPR1, ENVE5:REPR1, ENVE6:REPR1, ENVE7:REPR1, ENVE8:REPR1, ENVE9:REPR1, ENVE1:REPR10, ENVE2:REPR10, ENVE3:REPR10, ENVE5:REPR10, ENVE6:REPR10, ENVE7:REPR10, ENVE8:REPR10, ENVE9:REPR10, ENVE1:REPR11, ENVE2:REPR11, ENVE3:REPR11, ENVE5:REPR11, ENVE6:REPR11, ENVE7:REPR11, ENVE8:REPR11, ENVE9:REPR11, ENVE1:REPR12, ENVE2:REPR12, ENVE3:REPR12, ENVE5:REPR12, ENVE6:REPR12, ENVE7:REPR12, ENVE8:REPR12, ENVE9:REPR12, ENVE1:REPR13, ENVE2:REPR13, ENVE3:REPR13, ENVE4:REPR13, ENVE6:REPR13, ENVE7:REPR13, ENVE8:REPR13, ENVE9:REPR13, ENVE1:REPR14, ENVE2:REPR14, ENVE3:REPR14, ENVE4:REPR14, ENVE6:REPR14, ENVE7:REPR14, ENVE8:REPR14, ENVE9:REPR14, ENVE1:REPR15, ENVE2:REPR15, ENVE3:REPR15, ENVE4:REPR15, ENVE6:REPR15, ENVE7:REPR15, ENVE8:REPR15, ENVE9:REPR15, ENVE1:REPR16, ENVE2:REPR16, ENVE3:REPR16, ENVE4:REPR16, ENVE5:REPR16, ENVE7:REPR16, ENVE8:REPR16, ENVE9:REPR16, ENVE1:REPR17, ENVE2:REPR17, ENVE3:REPR17, ENVE4:REPR17, ENVE5:REPR17, ENVE7:REPR17, ENVE8:REPR17, ENVE9:REPR17, ENVE1:REPR18, ENVE2:REPR18, ENVE3:REPR18, ENVE4:REPR18, ENVE5:REPR18, ENVE7:REPR18, ENVE8:REPR18, ENVE9:REPR18, ENVE1:REPR19, ENVE2:REPR19, ENVE3:REPR19, ENVE4:REPR19, ENVE5:REPR19, ENVE6:REPR19, ENVE8:REPR19, ENVE9:REPR19, ENVE2:REPR2, ENVE3:REPR2, ENVE4:REPR2, ENVE5:REPR2, ENVE6:REPR2, ENVE7:REPR2, ENVE8:REPR2, ENVE9:REPR2, ENVE1:REPR20, ENVE2:REPR20, ENVE3:REPR20, ENVE4:REPR20, ENVE5:REPR20, ENVE6:REPR20, ENVE8:REPR20, ENVE9:REPR20, ENVE1:REPR21, ENVE2:REPR21, ENVE3:REPR21, ENVE4:REPR21, ENVE5:REPR21, ENVE6:REPR21, ENVE8:REPR21, ENVE9:REPR21, ENVE1:REPR22, ENVE2:REPR22, ENVE3:REPR22, ENVE4:REPR22, ENVE5:REPR22, ENVE6:REPR22, ENVE7:REPR22, ENVE9:REPR22, ENVE1:REPR23, ENVE2:REPR23, ENVE3:REPR23, ENVE4:REPR23, ENVE5:REPR23, ENVE6:REPR23, ENVE7:REPR23, ENVE9:REPR23, ENVE1:REPR24, ENVE2:REPR24, ENVE3:REPR24, ENVE4:REPR24, ENVE5:REPR24, ENVE6:REPR24, ENVE7:REPR24, ENVE9:REPR24, ENVE1:REPR25, ENVE2:REPR25, ENVE3:REPR25, ENVE4:REPR25, ENVE5:REPR25, ENVE6:REPR25, ENVE7:REPR25, ENVE8:REPR25, ENVE1:REPR26, ENVE2:REPR26, ENVE3:REPR26, ENVE4:REPR26, ENVE5:REPR26, ENVE6:REPR26, ENVE7:REPR26, ENVE8:REPR26, ENVE1:REPR27, ENVE2:REPR27, ENVE3:REPR27, ENVE4:REPR27, ENVE5:REPR27, ENVE6:REPR27, ENVE7:REPR27, ENVE8:REPR27, ENVE2:REPR3, ENVE3:REPR3, ENVE4:REPR3, ENVE5:REPR3, ENVE6:REPR3, ENVE7:REPR3, ENVE8:REPR3, ENVE9:REPR3, ENVE1:REPR4, ENVE3:REPR4, ENVE4:REPR4, ENVE5:REPR4, ENVE6:REPR4, ENVE7:REPR4, ENVE8:REPR4, ENVE9:REPR4, ENVE1:REPR5, ENVE3:REPR5, ENVE4:REPR5, ENVE5:REPR5, ENVE6:REPR5, ENVE7:REPR5, ENVE8:REPR5, ENVE9:REPR5, ENVE1:REPR6, ENVE3:REPR6, ENVE4:REPR6, ENVE5:REPR6, ENVE6:REPR6, ENVE7:REPR6, ENVE8:REPR6, ENVE9:REPR6, ENVE1:REPR7, ENVE2:REPR7, ENVE4:REPR7, ENVE5:REPR7, ENVE6:REPR7, ENVE7:REPR7, ENVE8:REPR7, ENVE9:REPR7, ENVE1:REPR8, ENVE2:REPR8, ENVE4:REPR8, ENVE5:REPR8, ENVE6:REPR8, ENVE7:REPR8, ENVE8:REPR8, ENVE9:REPR8, ENVE1:REPR9, ENVE2:REPR9, ENVE4:REPR9, ENVE5:REPR9, ENVE6:REPR9, ENVE7:REPR9, ENVE8:REPR9, ENVE9:REPR9.
## Interpret type III hypotheses with care.## Evaluating trait yield |=========================================| 100% 00:00:04 ## Method: REML/BLUP## Random effects: GEN, GEN:ENV## Fixed effects: ENV, REP(ENV)## Denominador DF: Satterthwaite's method## ---------------------------------------------------------------------------
## P-values for Likelihood Ratio Test of the analyzed traits
## ---------------------------------------------------------------------------
## model yield
## COMPLETE NA
## GEN 4.21e-22
## GEN:ENV 3.87e-15
## ---------------------------------------------------------------------------
## All variables with significant (p < 0.05) genotype-vs-environment interaction# Plotting BLUP
plot_blup(BLUP, which = "ge")