Workshop LME Bagian 2

# Linear Mixed Models (Multi-Level Modeling)
## Menggunakan `jamovi` (GAMLj): Bagian 2
 
### Rizqy Amelia Zein
* Dosen, [Fakultas Psikologi, Universitas Airlangga](https://psikologi.unair.ac.id)
* Anggota, [#SainsTerbuka Airlangga](https://sainsterbukaua.github.io/) <svg style="height:0.8em;top:.04em;position:relative;" viewBox="0 0 512 512"><path d="M256 48c114.953 0 208 93.029 208 208 0 114.953-93.029 208-208 208-114.953 0-208-93.029-208-208 0-114.953 93.029-208 208-208m0-40C119.033 8 8 119.033 8 256s111.033 248 248 248 248-111.033 248-248S392.967 8 256 8zm0 56C149.961 64 64 149.961 64 256s85.961 192 192 192 192-85.961 192-192S362.039 64 256 64zm0 44c19.882 0 36 16.118 36 36s-16.118 36-36 36-36-16.118-36-36 16.118-36 36-36zm117.741 98.023c-28.712 6.779-55.511 12.748-82.14 15.807.851 101.023 12.306 123.052 25.037 155.621 3.617 9.26-.957 19.698-10.217 23.315-9.261 3.617-19.699-.957-23.316-10.217-8.705-22.308-17.086-40.636-22.261-78.549h-9.686c-5.167 37.851-13.534 56.208-22.262 78.549-3.615 9.255-14.05 13.836-23.315 10.217-9.26-3.617-13.834-14.056-10.217-23.315 12.713-32.541 24.185-54.541 25.037-155.621-26.629-3.058-53.428-9.027-82.141-15.807-8.6-2.031-13.926-10.648-11.895-19.249s10.647-13.926 19.249-11.895c96.686 22.829 124.283 22.783 220.775 0 8.599-2.03 17.218 3.294 19.249 11.895 2.029 8.601-3.297 17.219-11.897 19.249z"/></svg>
* Relawan, [INA-Rxiv](https://inarxiv.id) 
* Researcher-in-training, [Institute for Globally Distributed Open Research and Education (IGDORE)](https://igdore.org/)

---

# Quiz time! 😄
## [Klik disini untuk menuju laman kuis.](https://menti.com)

---

# Coba kita lihat lebih dekat...

* Coba kita lakukan inspeksi visual sekali lagi dataset-sekolah.omv

* Buat *scatterplot* dimana **mandiri** menjadi **Y-Axis**, sedangkan **neu**, **hi**, **trust** sebagai **X-Axis**

* Kemudian masukkan **idsekolah** pada kolom **Group**
  - Fungsinya, kita akan mendapatkan garis regresi untuk masing-masing sekolah

* Apa yang terjadi?

---

---

# Ternyata *covariance*nya bervariasi di tiap kelompok! 😱

* Yang menarik, tidak hanya *intercept*, *slope* pendapatan personal juga bervariasi di setiap sekolah
  - Berdasarkan analisis yang kita lakukan pagi tadi, disimpulkan bahwa **pendapatan keluarga dan kemandirian anak korelasinya negatif** 
  - Tapi bisakah **kesimpulan yang sama** kita tarik untuk Sekolah B dan D?

* Hati-hati *ecological fallacy*!
  - Terjadi ketika kita salah menyimpulkan suatu gejala yang skalanya individual, padahal yang dianalisis oleh peneliti sesusungguhnya fenomena yang berlaku pada skala yang lebih besar (kelompok atau sub-kelompok)
]

]

---

# Struktur sampel bersarang/berjenjang

---

# Apa yang harus dilakukan?

* Pura-pura *ga tahu*, kita abaikan saja dan langsung menggunakan regresi OLS, dengan atau tanpa informasi mengenai pengelompokan data sebagai variabel kontrol.

- Masalahnya, data/observasi kita sangat bergantung pada pengelompokan unit analisis
  - **Nah, lalu melanggar asumsi OLS dong?** (data/observasi dan residual harus independen)
  - Efeknya, *standard error* yang diestimasi oleh model terlalu kecil (karena mengabaikan varians dependen variabel yang ditentukan oleh kelompok)
  - Varians variabel dependen yang tidak bisa dijelaskan (residual) akan makin besar
  - Kesimpulan/inferensi yang ditarik menjadi tidak tepat, sehingga risiko terjadinya *type I error* meningkat.

* Gimana kalau pengelompokan (*group status*) dimasukkan aja dalam regresi OLS sebagai variabel moderator
  - Dengan begitu, estimasi *standard error* disesuaikan dengan menggunakan *marginal model*

- Estimasi *standard error* akan lebih presisi, **tetapi** kita tetap tidak bisa mengestimasi *between-group variance*
  
--

* Kalo di agregat? Jadi unit analisis yang tadinya individual, menjadi kelompok.

- Ukuran sampel menjadi lebih sedikit, sehingga *statistical power* menjadi lebih rendah❗

---

# *Fixed* dan *random effects*

## Model *fixed effects*

]

## Model *random effects*

]

---

# *Full model*

---

# Kovarians (korelasi) antara *random intercept* dan *random slopes* (σ*U0**U1*)

* **Nilainya positif**, maka semakin tinggi *intercept* akan diasosiasikan dengan kemiringan garis yang lebih curam/*slopes* yang lebih besar

* Misalnya, di sekolah yang **rata-rata pendapatan** keluarga inti perbulan siswanya **tinggi**, maka **korelasi** antara pendapatan per bulan dengan tingkat kemandirian siswa **akan menguat**.

* **Nilainya negatif**, maka semakin tinggi *intercept* akan diasosiasikan dengan kemiringan garis yang lebih landai/*slopes* yang lebih kecil

* Misalnya, di sekolah yang **rata-rata pendapatan** keluarga inti perbulan siswanya **tinggi**, maka **korelasi** antara pendapatan per bulan dengan tingkat kemandirian siswa **akan melemah**.
]
  
.pull-right[

]

---

# Parameter yang diestimasi dalam *linear mixed effects*

* *Fixed intercept* (*c**00*)

* *Fixed slopes* (*c**10*)

* Varians *random intercept* (σ2*U0*)

* Varians *random slopes* (σ2*U1*)

* Kovarians (korelasi) antara *random intercept* dan *random slopes* (σ*U0**U1*)

* Varians residual level-1 (σ2*e*)

---

# Yuk kita coba! ✊
## Pastikan *module* GAMLj sudah terpasang di jamovi

![](https://media.giphy.com/media/cjudIH6tZrvcPD2P4j/giphy.gif)

---

# Latihan (3): Kembali ke sekolah 🏫

Setelah menginspeksi data secara visual, kita tahu bahwa korelasi antara **pendapatan keluarga** dengan **tingkat kemandirian anak** adalah yang paling bervariasi (daripada prediktor yang lainnya). Oleh karena itu, kita akan membuat *linear mixed model* dengan **pendapatan keluarga** sebagai prediktor, dan **tingkat kemandirian anak** sebagai variabel dependen.

## Bikin "model kosong"

* Yaitu model yang isinya hanya *intercept* saja, tidak ada prediktornya (*slopes*)

* Pada *menu bar*, klik **Linear Models**, pilih **mixed models**
  - Masukkan **mandiri** dalam kolom **dependent variable**
  - Masukkan **idsekolah** dalam kolom **cluster variables**
  - Pada menu **random effects** masukkan **intercept|idsekolah** dalam kolom **random coefficients**
  
* Catat nilai AIC yang tersedia dalam tabel **model info**

* Yuk ingat-ingat lagi, AIC fungsinya untuk apa ya? 😕

---

# Latihan (3)

## Bikin model dengan prediktor (*linear mixed model*)

* Masukkan **hi** dalam kolom **covariates**

* Pada menu **random effects**, masukkan juga **hi|idsekolah**, karena kita akan mengestimasi **random slopes**nya juga
  - Centang opsi **LRT for Random Test**

* Pada menu **covariates scaling**, ubah **centered** menjadi **cluster-based centered**
  - Berkaitan dengan **partitioning** (akan dijelaskan di bagian selanjutnya)

---

# *Fixed coefficients*

* Tes kecocokan model (*Omnibus Test*) tidak signifikan (*F*(1,3.98)=0.309, p=.608) 
  - Berbeda sekali dengan hasil OLS
  - **Kemungkinan besar** korelasi antara pendapatan keluarga per bulan dengan kemandirian **tidak linier**

* Kita tidak punya cukup bukti untuk menolak hipotesis nol, bahwa pendapatan keluarga per bulan dengan kemandirian tidak berkaitan (*B*=-0.207 95% CI [-0.937, 0.523], *SE*=0.372, *t*=-0.556, *p*=.608).
  - Sekali lagi, berbeda sekali dengan hasil OLS

]

]

---

# *Random coefficients*

* Varians kemandirian paling banyak dijelaskan oleh rata-rata tingkat kemandirian anak di masing-masing sekolah (varians *random intercept* (σ2*U0*)) daripada oleh varians *random slopes* (σ2*U1*).

]

]

---

# *Random coefficients*

* Menguji efek sekolah (kelompok)
  - *Intra-class correlation*, yaitu merupakan proporsi total varians variabel dependen yang dapat dijelaskan oleh variasi antar kelompok
  
  - *Likelihood ratio test* (LRT), yaitu teknik untuk menguji ada/tidaknya perbedaan varians antar-kelompok
  
  - Keduanya juga bisa berfungsi sebagai indikator perlu/tidaknya `lme` dilakukan
  
]

]

---

# *Random coefficients*

* Menguji efek sekolah (kelompok) / (*testing group effects*)

- ICC=0.931, artinya 93.1% varians tingkat kemandirian anak dijelaskan oleh perbedaan kelompok. ICC diatas 0.1 biasanya menunjukkan `lme` adalah opsi yang lebih baik daripada OLS.
 
 - LRT juga menunjukkan bahwa kita dapat menolak hipotesis bahwa tidak adanya perbedaan varians tingkat kemandirian antar-kelompok (LRT(2)=74.0, *p*<.001).
]

]

---

# *Model Comparison*

* AIC
  - Apabila kita membandingkan "model kosong" dengan model yang ada prediktor, maka model yang terakhir lebih mampu menjelaskan varians kemandirian anak.

* R2 ([Nakagawa & Schielzeth, 2012](https://besjournals.onlinelibrary.wiley.com/doi/full/10.1111/j.2041-210x.2012.00261.x))
 
 - *Marginal*, merupakan proporsi varians variabel dependen yang dapat dijelaskan oleh **hanya fixed models** saja
 
 - *Conditional*, merupakan proporsi varians variabel dependen yang dapat dijelaskan oleh **fixed** dan **random models** sekaligus
 
 - Varians tingkat kemandirian anak yang dapat dijelaskan oleh *fixed model* saja hanya 0.12%, sedangkan yang dapat dijelaskan oleh keseluruhan model adalah 93.18%.

]

]

---

# *Contextual effects* & *partitioning*

* *Within-group effect*
  - Seberapa besar selisih Y dari 2 orang yang berada di **kelompok yang sama**, ketika **selisih Xnya** sebesar 1 poin?
  - Seberapa besar perbedaan tingkat kemandirian dua orang anak yang berada dalam **sekolah yang sama**, ketika selisih **tingkat pendapatan keluarga** mereka berbeda sebesar 1 poin?

* *Between-group effect*
  - Seberapa besar selisih Y dari 2 orang yang berbeda (dari 2 kelompok yang berbeda), namun berada dalam **posisi relatif yang sama** (dibandingkan dengan rerata kelompok), ketika **selisih X** mereka sebesar 1 poin?
  - Seberapa besar perbedaan tingkat kemandirian 2 orang anak yang berasal dari **sekolah yang berbeda**, namun pada **posisi yang sama** apabila **dibandingkan dengan rata-rata** kemandirian anak-anak di **sekolah mereka masing-masing** (misal sama-sama yang paling/paling tidak manja), apabila **selisih pendapatan keluarga** mereka = 1 poin?

* *Contextual effect*
  - Seberapa besar selisih Y dua orang dari kelompok yang berbeda, namun dengan **X yang sama**, ketika **rerata X kelompoknya** berbeda sebesar 1 poin.
  - Seberapa besar selisih tingkat kemandirian 2 orang anak dari **sekolah yang berbeda**, namun memiliki **tingkat pendapatan keluarga yang sama** persis, ketika **rata-rata tingkat pendapatan keluarga anak-anak di sekolahnya** berbeda sebesar 1 poin?

---

# *Contextual effects* & *partitioning*

* Untuk menghitung *contextual effect*, kita harus melakukan *partitioning* terlebih dahulu

* Umumnya yang dipartisi/*centering* adalah **variabel X, bukan Y**

* *Group-mean centering*
  - Nilai X individu dikurangi rata-rata X kelompoknya
  - Pendapatan keluarga anak A dikurangi rata-rata pendapatan keluarga anak-anak di sekolahnya

* *Grand-mean centering*
  - Nilai X individu dikurangi rata-rata X pada seluruh sampel
  - Pendapatan keluarga anak A dikurangi rata-rata pendapatan keluarga seluruh anak yang menjadi sampel

* *Contextual effect* = *Between-group effect* - *Within-group effect*
  - *Contextual effect* yang positif artinya **kelompok** dengan **rata-rata X yang lebih tinggi**, cenderung memiliki *intercept* (rata-rata Y) yang **lebih tinggi**, atau begitu pula sebaliknya.
  - *Contextual effect* yang negatif artinya **kelompok** dengan **rata-rata X yang lebih tinggi**, cenderung memiliki *intercept* (rata-rata Y) yang **lebih rendah**, atau begitu pula sebaliknya.

---

# Latihan (4): *contextual effects*

* Lakukan `lme` dengan memasukkan **hi_group_centered** dan **hi_grand_mean_centered** dalam satu model yang sama

* Masukkan kedua variabel tersebut dalam **fixed coefficients** dan **random coefficients**

* Lihat *fixed slopes*-nya untuk kedua prediktor

]

]

---

# *Contextual effects*

* Kita tidak punya bukti yang meyakinkan bahwa kita dapat menolak hipotesis bahwa *within* (*B*=0.257 95% CI [-0.115, 0.629], *SE*=0.190, *t*=1.354, *p*=.213), maupun *between-group effect* (*B*=-0.454 95% CI [-1.207, 0.299], *SE*=0.384, *t*=-1.182, *p*=.315) tidak dapat menjelaskan varians tingkat kemandirian anak.

* *Contextual effects* = -0.711

]

]

---

# Bagaimana melaporkannya? (1)

"...untuk menguji hipotesis bahwa ada perbedaan rerata tingkat kemandirian anak, dan korelasi antara pendapatan keluarga dengan tingkat kemandirian anak di masing-masing sekolah, peneliti melakukan analisis *linear mixed effect*.

Tingkat kemandirian anak dijelaskan sebagai fungsi dari tingkat pendapatan keluarga, dengan mengontrol asal sekolah (PAUD) anak. Sebelum melakukan analisis, tingkat pendapatan keluarga dipartisi dengan cara menguranginya dengan rata-rata tingkat pendapatan keluarga di masing-masing sekolah (*group-mean/cluster-based centering*).

Pengujian model menghasilkan kesimpulan bahwa model tidak cocok menggambarkan data (*F*(1,3.98)=0.309, p=.608), sehingga kami gagal menolak hipotesis bahwa tingkat pendapatan keluarga dan kemandirian anak, tidak berkorelasi. Namun, ada kemungkinan korelasi antara tingkat pendapatan keluarga dan kemandirian anak tidak linier.

Model *fixed effects* menunjukkan bahwa peneliti tidak punya cukup bukti untuk menolak hipotesis nol, bahwa pendapatan keluarga per bulan dengan kemandirian tidak berkaitan (*B*=-0.207 95% CI [-0.937, 0.523], *SE*=0.372, *t*=-0.556, *p*=.608).

---

# Bagaimana melaporkannya? (2)

"...model *random effects* menunjukkan bahwa peneliti dapat menolak hipotesis bahwa tidak ada perbedaan varians tingkat kemandirian antar-kelompok (LRT(2)=74.0, *p*<.001). Varians kemandirian anak paling banyak dijelaskan oleh rata-rata tingkat kemandirian anak di masing-masing sekolah (varians *random intercept* (σ2*U0*)) daripada oleh varians *random slopes* (σ2*U1*). Selain itu, 93.1% varians tingkat kemandirian anak dijelaskan oleh perbedaan kelompok (ICC=0.931).

Peneliti tidak memiliki cukup bukti yang meyakinkan untuk menolak hipotesis yang menyatakan bahwa *within* (*B*=0.257 95% CI [-0.115, 0.629], *SE*=0.190, *t*=1.354, *p*=.213), maupun *between-group effect* (*B*=-0.454 95% CI [-1.207, 0.299], *SE*=0.384, *t*=-1.182, *p*=.315) tidak dapat menjelaskan varians tingkat kemandirian anak.

*Contextual effects* ditemukan sebesar -0.711, artinya, anak yang bersekolah di dua tempat yang berbeda, dengan selisih rata-rata tingkat kemandirian anak-anak di dua sekolah tersebut sebesar 1 poin, maka tingkat kemandirian mereka berbeda sebesar -0.711 poin, apabila diasumsikan keluarga mereka memiliki pendapatan yang sama besarnya. Tanda negatif dari *contextual effect* mengindikasikan bahwa sekolah yang rata-rata pendapatan keluarga siswanya lebih tinggi, cenderung memiliki rata-rata tingkat kemandirian yang rendah..."

---

# Latihan mandiri (2)

* Lakukan analisis `lme` untuk mengetahui:

- Apakah varians **tingkat kemandirian anak** dapat dijelaskan oleh sekolah tempat anak tersebut belajar?
  
  - Apakah varians korelasi antara kecenderungan **neuroticism** ibu dengan **kemandirian anak** juga dapat dijelaskan oleh sekolah tempat anak tersebut belajar?
  
  - Seberapa besar perbedaan **tingkat kemandirian** dua orang anak yang berada di **sekolah yang berbeda**, yang **ibunya sama-sama pencemas**, apabila **rata-rata kecemasan** wali murid di **dua sekolah tersebut** berbeda sebesar 1 poin?

## [Klik disini untuk mengakses lembar kerja](https://forms.gle/87BrLEdC321khd8g6)

---

# Yang belum dibahas...

* Kalau korelasi antara X dan Y tidak linier, pakai apa dong?
  - Jelas tidak bisa menggunakan `lme`. Alternatifnya, bisa menggunakan [*generalised additive model* (GAM)](https://en.wikipedia.org/wiki/Generalized_additive_model).

* Kalau prediktornya level-2, bagaimana?

* Bisa ga `lme` digunakan untuk mengestimasi perubahan Y pada *time-series*?

* Gimana cara merencanakan jumlah sampelnya?

* Gimana kalo sampelnya bersarang/berjenjang level-3, bahkan lebih?

* Gimana kalo terjadi interaksi antara variabel prediktor level-1 dengan level-2 (*cross-level interactions*)?

]

]

---

# The problem with linear relationship

---

# Terima kasih banyak! 😉

![](https://media.giphy.com/media/hrBSJ2So6iTo4/giphy.gif)

Paparan disusun dengan menggunakan <svg style="height:0.8em;top:.04em;position:relative;" viewBox="0 0 581 512"><path d="M581 226.6C581 119.1 450.9 32 290.5 32S0 119.1 0 226.6C0 322.4 103.3 402 239.4 418.1V480h99.1v-61.5c24.3-2.7 47.6-7.4 69.4-13.9L448 480h112l-67.4-113.7c54.5-35.4 88.4-84.9 88.4-139.7zm-466.8 14.5c0-73.5 98.9-133 220.8-133s211.9 40.7 211.9 133c0 50.1-26.5 85-70.3 106.4-2.4-1.6-4.7-2.9-6.4-3.7-10.2-5.2-27.8-10.5-27.8-10.5s86.6-6.4 86.6-92.7-90.6-87.9-90.6-87.9h-199V361c-74.1-21.5-125.2-67.1-125.2-119.9zm225.1 38.3v-55.6c57.8 0 87.8-6.8 87.8 27.3 0 36.5-38.2 28.3-87.8 28.3zm-.9 72.5H365c10.8 0 18.9 11.7 24 19.2-16.1 1.9-33 2.8-50.6 2.9v-22.1z"/></svg> *package* [**xaringan**](https://github.com/yihui/xaringan) dengan *template* dan *fonts* dari `R-Ladies`.

*Chakra* dibuat dengan [remark.js](https://remarkjs.com), [**knitr**](http://yihui.name/knitr), dan [R Markdown](https://rmarkdown.rstudio.com).