Model Campuran Linear Terampat/ Generalized Linear Mixed Model (GLMM)

Author

Achmad Fauzan

1 Pengantar

Pada Model Linear (Linear Model, LM), seperti regresi linear, memerlukan asumsi bahwa peubah respon \(Y\) menyebar normal.
Pada kenyataannya banyak ditemukan bahwa peubah respon \(Y\) tidak menyebar Normal. Misalnya menyebar Binomial, Binomial-negatif, Multinomial, Poisson, Gammar, Eksponensial, Chi-Square, Log-Normal, Beta, dsb.
Dikembangkan Model Linear Terampat (Generalized Linear Model, GLM) untuk mengatasi masalah ini.
Kemudian GLM ini diperluas menjadi Model Campuran Linear Terampat (Generalized Linear Mixed Model, GLMM) dengan menambahkan pengaruh acak (random effect).

Kalau diilustrasikan, dapat dilihat pada Gambar Figure 1

Linear Model (LM) >> Model regresi linear klasik tanpa efek acak.
Generalized Linear Model (GLM) >> Perluasan LM untuk menangani distribusi non-normal (keluarga eksponensial, misalnya binomial, poisson, dsb) pada variabel responnya.
Linear Mixed Model (LMM) >> Perluasan LM dengan efek acak (random effects) untuk menangani data hierarkis atau berulang.
Generalized Linear Mixed Model (GLMM) >> Kombinasi GLM dan LMM, yang menangani distribusi non-normal dan efek acak.

2 Linear Mixed Models(LMM)

Mixed-Effect disebut “Mixed” karena menggabungkan fixed effect dan random effect secara bersamaan.
Tidak ada aturan secara objektif kapan efek itu masuk efek acak atau efek tetap, karena apa yang masuk sebagai efek acak disuatu penelitian bisa jadi menjadi efek tetap dipenelitian yang lain, sehingga bergantung penelitiannya. Namun ada beberapa pedoman, lebih kurangnya sebagai berikut.
Fixed effect merepresentasikan efek pada Tingkat populasi (rata-rata efek yang berlaku secara umum dan konsisten dalam berbagai eksperimen). Biasanya fixed effect menjadi parameter yang ingin kita estimasi atau kita duga duga dalam regresi linier.
Random effect menggambarkan seberapa besar variasi tren ini diantara kelompok (misalnya, peserta atau item dalam penelitian). Bisa juga dikatakan Random effect adalah parameter yang bisa berbeda antar kelompok data yang saling berkaitan. Jika mencoba mengontrol ketidak mandirian yang dihasilkan oleh struktur data bersarang, dimana kita tidak benar-benar peduli pada faktornya sendiri tetapi ingin memperhitungkan perbedaan potensial antara kelompok akan menjadi efek acak (dan ini merupakan situasi paling umum pada efek campuran). Misalnya, jika kita mengukur sesuatu beberapa kali pada orang yang sama, rata-rata dari hasil pengukuran itu bisa menjadi satu perkiraan nilai. Setiap individu bisa memiliki perkiraan nilai yang berbeda.
Dengan kata lain, fixed effects menangkap pola rata-rata, sementara random effects menangkap variasi antar kelompok dalam data.

Sebagai ilustrasi, berikut 3 contoh yang relevan berkaitan dengan fixed effects dan random effect

Studi tentang prestasi siswa di sekolah

Seorang peneliti ingin mengetahui pengaruh jumlah jam belajar terhadap nilai ujian siswa di berbagai sekolah.

Fixed Effect: Jumlah jam belajar (karena efeknya diharapkan sama untuk semua siswa).
Random Effect: Sekolah (karena tiap sekolah mungkin memiliki standar pengajaran yang berbeda).
Model ini mengukur hubungan antara jam belajar dan nilai ujian, tetapi juga mempertimbangkan bahwa setiap sekolah mungkin memiliki pengaruh berbeda terhadap nilai siswa.

Pengaruh obat terhadap tekanan darah

Sebuah studi klinis meneliti efektivitas obat terhadap tekanan darah pasien dari berbagai rumah sakit.

Fixed Effect: Jenis obat yang diberikan.
Random Effect: Rumah sakit (karena setiap rumah sakit bisa memiliki faktor lingkungan atau prosedur medis yang berbeda).
Model ini mengevaluasi apakah obat memiliki efek signifikan terhadap tekanan darah, sekaligus menangkap perbedaan antar rumah sakit.

Studi tentang kecepatan lari atlet

Seorang pelatih ingin mengetahui hubungan antara pola latihan dan kecepatan lari atlet di beberapa klub olahraga.

Fixed Effect: Pola latihan (karena efek latihan dianggap berlaku untuk semua atlet).
Random Effect: Klub olahraga (karena setiap klub bisa memiliki fasilitas atau metode pelatihan yang berbeda).
Model ini mengukur pengaruh pola latihan terhadap kecepatan lari, sambil mempertimbangkan variasi antar klub olahraga.

Oke, kita lanjutkan ya ke contoh.

Misalnya kita memiliki data berat badang (Kg) 4 orang mencoba diet tertentu. Kolom kedua menyatakan berat orang tersebut sebelum diet, sementara kolom ketiga dan keempat menyatakan berat badan setelah melakukan diet. Kita ingin mengetahui Penurunan berat badan per minggu.

Person	Diet	Before the diet	After 1 week	After 2 weeks	After 3 weeks
Person 1	A	102	97	95	93
Person 2	A	96	93	87	85
Person 3	B	83	79	78	74
Person 4	B	79	77	75	72

Sumber: Linear mixed effects models

Apabila divisualisakan data diatas dapat disajikan pada gambar Figure 2

Pada Figure 2, warna pada titik menyatakan individu, sumbu horizontal (\(X\)) menyatakan kurun waktu (minggu) dan sumbu \(Y\) menyatakan berat badan (Kg). Apabila kita jalankan menggunakan program R, berikut langkahnya.

Data yang digunakan dapat dituliskan sebagai berikut.

rm(list=ls())   # Menghapus semua dataset dan variabel
graphics.off()  # Menutup semua gambar

#Data
df <- data.frame(
  Weight = c(102,97, 95, 93, 96, 93, 87, 85, 83, 79, 78, 74, 79, 77, 75, 72),
  Subjects = rep(1:4, each = 4),
  Diet = rep(c("A", "B"), each = 8),
  Weeks = rep(c(0, 1, 2, 3), times = 4)
)

# Tampilkan data frame
print(df)

   Weight Subjects Diet Weeks
1     102        1    A     0
2      97        1    A     1
3      95        1    A     2
4      93        1    A     3
5      96        2    A     0
6      93        2    A     1
7      87        2    A     2
8      85        2    A     3
9      83        3    B     0
10     79        3    B     1
11     78        3    B     2
12     74        3    B     3
13     79        4    B     0
14     77        4    B     1
15     75        4    B     2
16     72        4    B     3

Misalkan dari data tersebut kita lakukan regresi linear langsung (tanpa memperhatikan subjects) dengan variabel prediktor adalah Weeks dan variabel respon adalah Weight.

# Model regresi linear
model <- lm(Weight ~ Weeks , data = df)

# Ringkasan hasil regresi
summary(model)


Call:
lm(formula = Weight ~ Weeks, data = df)

Residuals:
    Min      1Q  Median      3Q     Max 
-10.775  -8.056  -1.325   7.219  12.225 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   89.775      3.851  23.312 1.33e-12 ***
Weeks         -2.975      2.058  -1.445     0.17    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 9.206 on 14 degrees of freedom
Multiple R-squared:  0.1298,    Adjusted R-squared:  0.06767 
F-statistic: 2.089 on 1 and 14 DF,  p-value: 0.1704

Model regresi linear sederhana diatas, dapat dituliskan menjadi Persamaan \(\eqref{eq:regresi}\)

\[ \begin{equation} \begin{aligned} \hat{Y} &= \beta_0 + \beta_1 \cdot X \\ &= 89.775 -2.975 \cdot X \end{aligned} \tag{1} \label{eq:regresi} \end{equation} \]

\(Y\): berat badan dan \(X\): minggu.

Persamaan \(\eqref{eq:regresi}\) menunjukkan bahwa rata-rata setiap minggu, berat badan mengalami penurunan sebesar 2.975 kg. Namun jika kita lihat nilai \(p-val = 0.17 > 0.05 = \alpha\) artinya kemiringan tidak berbeda signifikan dari nol (\(H_0: slope (\beta_0) =0\)).

Karena kita memiliki pengukuran berulang pada individu yang sama, model efek linier campuran lebih tepat.

2.1 LMM (intecept acak/beda, slope sama)

Pada contoh pertama, kita akan menganalisis model dengan intercept acak yang bearti kita akan memperkirakan intercept untuk setiap individu, tetapi kita mengasumsikan semua individu memiliki kemiringan (slope) yang sama.

Notasinya ditulis pada notasi \(\eqref{eq:lme1}\) (ilustrasi dari package lme4).

\[ \begin{equation} \begin{aligned} \text{Weight} \sim \text{Weeks} + (1|\text{Subjects}) \end{aligned} \tag{2} \label{eq:lme1} \end{equation} \]

Pada notasi \(\eqref{eq:lme1}\) , efek tetap diletakan pada Weeks, sementara efek acak adalah Subjects.

Disebelah kiri tanda \(|\), digunakan untuk menggunakan random intercept atau slope, atau keduanya.

1: efek acak intersep (Random intercept only), Subjek memiliki intersep yang berbeda-beda, tetapi slope (Weeks) tetap sama untuk semua subjek. Modelnya dituliskan pada notasi \(\eqref{eq:lme2}\)

\[ \begin{equation} \begin{aligned} \text{Y}_{ij} = \beta_0 + b_{0j} + \beta_1 X +\epsilon_{ij} \end{aligned} \tag{2} \label{eq:lme2} \end{equation} \]

\(i\): indeks pengamatan individu dalam suatu subjek (minggu), \(j\): kelompok/ individu yang berbeda

0: Efek acak slope (Random Slope Only), Subjek memiliki slope yang berbeda-beda, tetapi intersep semua subjek sama. Modelnya dituliskan pada notasi \(\eqref{eq:lme3}\)

\[ \begin{equation} \begin{aligned} \text{Y}_{ij} = \beta_0 + (\beta_1+b_{1j}) X +\epsilon_{ij} \end{aligned} \tag{3} \label{eq:lme3} \end{equation} \]

Efek Acak Intersep dan Slope (Random Intercept and Slope), Baik intersep maupun slope bersifat acak, setiap subjek memiliki regresi sendiri. Modelnya dituliskan pada notasi \(\eqref{eq:lme4}\)

\[ \begin{equation} \begin{aligned} \text{Y}_{ij} = (\beta_0 + b_{0j}) (\beta_1+b_{1j}) X +\epsilon_{ij} \end{aligned} \tag{4} \label{eq:lme4} \end{equation} \]

Pada contoh diatas, kalau hanya menggunakan regresi linier (sederhana), kemiringan tidak berpengaruh signifikan (karena titik data jauh dari garis slope)

Namun jika menggunakan LLM akan diperoleh \(p-val<0.001 (H_0 \text{ditolak})\) dan dapat disimpulkan kemiringan (Weeks) berpengaruh secara signifikan. Oleh karena itu digunakan LMM untuk mengatasi variasi berat badan individu tersebut.

Sintaks berikut mengasumsikan Subjek memiliki intersep yang berbeda-beda, tetapi slope tetap sama untuk semua subjek.

# install.packages("lmerTest")
library(lmerTest)
M4<-lmer(Weight~Weeks+(1|Subjects),REML=T,data=df)
# Subjek memiliki intersep yang berbeda-beda, tetapi slope (Weeks) tetap sama untuk semua subjek.
summary(M4)

Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: Weight ~ Weeks + (1 | Subjects)
   Data: df

REML criterion at convergence: 66.2

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-1.53483 -0.73032 -0.01984  0.67606  1.14606 

Random effects:
 Groups   Name        Variance Std.Dev.
 Subjects (Intercept) 97.359   9.867   
 Residual              1.294   1.138   
Number of obs: 16, groups:  Subjects, 4

Fixed effects:
            Estimate Std. Error      df t value Pr(>|t|)    
(Intercept)  89.7750     4.9564  3.0359   18.11 0.000342 ***
Weeks        -2.9750     0.2544 11.0000  -11.69 1.52e-07 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
      (Intr)
Weeks -0.077

# Tambahan
# M4<-lmer(Weight ~ Weeks + (0 + Weeks | Subjects),REML=F,data=df) 
# Subjek memiliki slope yang berbeda-beda, tetapi intersep semua subjek sama.

# M4<-lmer(Weight ~ Weeks + (1 + Weeks | Subjects),REML=F,data=df) 
# Baik intersep maupun slope bersifat acak, setiap subjek memiliki regresi sendiri.

Hasil diatas menampilkan estimasi dari fixed effect (estimasi keseluruhan/ rata-rata), sementara estimasi dari efek acaknya dapat dikeluarkan dari sintaks berikut.

ranef(M4)

$Subjects
  (Intercept)
1   11.399612
2    4.921144
3   -6.789933
4   -9.530824

with conditional variances for "Subjects"

Nilai dari efek acak mewakili seberapa besar intercept dari subject ke-j dibandingkan dengan intercept keseluruhan.

Persamaan modelnya disajikan pada Persamaan \(\eqref{eq:lme5}\)

\[ \begin{equation} \begin{aligned} \hat {Y_{ij}} = \beta_0 + b_{0j} + \beta_1 X \end{aligned} \tag{5} \label{eq:lme5} \end{equation} \]

\(i\): Minggu ke-\(i\) \(j\): Subject ke-\(j\)

Sebagai ilustrasi sebagai berikut.

Prediksi berat badan pada minggu ke 1, subject ke-2

\[ \begin{equation} \begin{aligned} \hat {Y_{12}} &= \beta_0 + b_{02} + \beta_1 X \\ &= 89.775 + 4.921144 - (1) 2.9750 \\ &= 91.721 \end{aligned} \end{equation} \]

intercept variabel acak dari subject ke-2 4.921144 berarti intersep pada Subject 2 lebih tinggi sebesar 4.921144 jika dibandingkan dengan intercept keseluruhan sehingga intercept pada subject ke-2 adalah \(89.775+ 4.921144 = 94.696\)

Prediksi berat badan pada minggu ke 2, Subject ke-3

\[ \begin{equation} \begin{aligned} \hat {Y_{23}} &= \beta_0 + b_{03} + \beta_1 X \\ &= 89.775 -6.789933 - 2 (2.9750) \\ &= 70.245 \end{aligned} \end{equation} \] intercept variabel acak dari subject ke-3 -6.789933 berarti intersep pada Subject 3 lebih rendah sebesar 6.789933 jika dibandingkan dengan intercept keseluruhan sehingga intercept pada subject ke-3 adalah \(89.775-6.789933 = 76.195\)

2.2 LMM (intecept acak/beda, slope juga acak/beda)

Selanjutnya, dieksplorasi kondisi jika intercept maupun slope dianggap acak, artinya setiap subjek memiliki regresi sendiri.

M5<-lmer(Weight ~ Weeks + (1 + Weeks | Subjects),REML=F,data=df) 
summary(M5)

Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
  method [lmerModLmerTest]
Formula: Weight ~ Weeks + (1 + Weeks | Subjects)
   Data: df

     AIC      BIC   logLik deviance df.resid 
    80.1     84.8    -34.1     68.1       10 

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-1.66980 -0.49134  0.01578  0.59857  1.34248 

Random effects:
 Groups   Name        Variance Std.Dev. Corr 
 Subjects (Intercept) 81.8806  9.0488        
          Weeks        0.1494  0.3865   -0.88
 Residual              0.9375  0.9683        
Number of obs: 16, groups:  Subjects, 4

Fixed effects:
            Estimate Std. Error      df t value Pr(>|t|)    
(Intercept)  89.7750     4.5425  4.0000   19.76 3.87e-05 ***
Weeks        -2.9750     0.2902  4.0000  -10.25  0.00051 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
      (Intr)
Weeks -0.634

ranef(M5)

$Subjects
  (Intercept)      Weeks
1   11.910177 -0.3533448
2    5.428958 -0.3185862
3   -7.159986  0.2485626
4  -10.179149  0.4233684

with conditional variances for "Subjects"

Terdapat korelasi negatif antara kemiringan random slopes and intercepts sebesar-0.88, atau dalam artian individu dengan berat awal yang lebih besar menurun lebih cepat jika dibandingkan individu yang memiliki berat awal yang lebih ringan.

2.3 LMM (intecept tetap, slope acak)

Ini hanya sebatas ekplorasi, ketika intercept dianggap tetap namun slope acak.

M6<-lmer(Weight ~ Weeks + (0 + Weeks | Subjects),REML=F,data=df)

3 Generalized Linear Mixed Models (GLMM)

Dalam GLMM, peubah respon diasumsikan mempunyai sebaran termasuk ke dalam keluarga eksponensial (exponential family)/ tidak hanya sebaran normal, yaitu:

\[ f_Y(y; \theta, \phi) = \exp \left\{ \frac{y\theta - b(\theta)}{a(\phi)} + c(y, \phi) \right\} \]

3.1 Contoh

Berikut adalah contoh GLMM dengan variabel \(Y\) memiliki sebaran binomial.

Data yang digunakan bersumber dari Penelitian dari Crowder (1978) menunjukkan hasil dari eksperimen faktorial 2×2 yang melibatkan dua faktor utama:

Varietas Benih (seed Variety) (\(X_1\))
Jenis ekstrak akar (Type of Root Extract) (\(X_2\))

Hasil yang diamati adalah jumlah benih yang berhasil berkecambah (germinated) dari sejumlah \(n\) benih yang ditanam. Dengan kata lain, setiap pengamatan mencatat jumlah benih yang berkecambah dalam satu kelompok, sementara \(n_i\) adalah total benih yang ditanam dalam kelompok tersebut.

Model yang umum digunakan untuk data seperti ini adalah distribusi binomial, yaitu:

\[Y_i \sim \text{Binomial} (n_i,p_i)\]

di mana \(p_i\) adalah probabilitas keberhasilan (perkecambahan) yang dapat dipengaruhi oleh faktor varietas benih dan jenis ekstrak akar.

Datanya adalah sebagai berikut DATA dan data tersebut berasal dari jurnal berikut JURNAL

\[\begin{array}{|c|c|c|c|c|} \hline \textbf{plate} & \textbf{seed} & \textbf{extract} & \textbf{r} & \textbf{n} \\ \hline 1 & 075 & \text{Bean} & 10 & 39 \\ 2 & 075 & \text{Bean} & 23 & 62 \\ 3 & 075 & \text{Bean} & 23 & 81 \\ 4 & 075 & \text{Bean} & 26 & 51 \\ 5 & 075 & \text{Bean} & 17 & 39 \\ 6 & 073 & \text{Bean} & 8 & 16 \\ 7 & 073 & \text{Bean} & 10 & 30 \\ 8 & 073 & \text{Bean} & 8 & 28 \\ 9 & 073 & \text{Bean} & 23 & 45 \\ 10 & 073 & \text{Bean} & 0 & 4 \\ 11 & 075 & \text{Cucumber} & 5 & 6 \\ 12 & 075 & \text{Cucumber} & 53 & 74 \\ 13 & 075 & \text{Cucumber} & 55 & 72 \\ 14 & 075 & \text{Cucumber} & 32 & 51 \\ 15 & 075 & \text{Cucumber} & 46 & 79 \\ 16 & 075 & \text{Cucumber} & 10 & 13 \\ 17 & 073 & \text{Cucumber} & 3 & 12 \\ 18 & 073 & \text{Cucumber} & 22 & 41 \\ 19 & 073 & \text{Cucumber} & 15 & 30 \\ 20 & 073 & \text{Cucumber} & 32 & 51 \\ 21 & 073 & \text{Cucumber} & 3 & 7 \\ \hline \end{array}\]

Data tersebut memiliki 5 kolom dan 21 baris. Penjelasannya adalah sebagai berikut.

plate: Plate number/ nomor pelat tempat percobaan dilakukan
seed: Jenis benih yang digunakan, dengan 2 kategori:
1. O75 (O. aegyptiaca 75)
2. O73 (O. aegyptica 73)
extract: Jenis ekstrak akar yang digunakan dalam percobaan, dengan 2 kategori:
1. Bean (ekstrak akar kacang)
2. Cucumber (ekstrak akar mentimun)
r: Jumlah benih yang berhasil berkecambah pada masing-masing plate
n: Total jumlah benih yang ditanam pada masing-masing plate

library(hglm.data)

# Memuat dataset "seeds"
data(seeds, package = "hglm.data")

# Menampilkan beberapa baris pertama dari dataset
seeds

   plate seed  extract  r  n
1      1  O75     Bean 10 39
2      2  O75     Bean 23 62
3      3  O75     Bean 23 81
4      4  O75     Bean 26 51
5      5  O75     Bean 17 39
6      6  O73     Bean  8 16
7      7  O73     Bean 10 30
8      8  O73     Bean  8 28
9      9  O73     Bean 23 45
10    10  O73     Bean  0  4
11    11  O75 Cucumber  5  6
12    12  O75 Cucumber 53 74
13    13  O75 Cucumber 55 72
14    14  O75 Cucumber 32 51
15    15  O75 Cucumber 46 79
16    16  O75 Cucumber 10 13
17    17  O73 Cucumber  3 12
18    18  O73 Cucumber 22 41
19    19  O73 Cucumber 15 30
20    20  O73 Cucumber 32 51
21    21  O73 Cucumber  3  7

summary(seeds)

     plate     seed        extract         r               n        
 Min.   : 1   O73:10   Bean    :10   Min.   : 0.00   Min.   : 4.00  
 1st Qu.: 6   O75:11   Cucumber:11   1st Qu.: 8.00   1st Qu.:16.00  
 Median :11                          Median :17.00   Median :39.00  
 Mean   :11                          Mean   :20.19   Mean   :39.57  
 3rd Qu.:16                          3rd Qu.:26.00   3rd Qu.:51.00  
 Max.   :21                          Max.   :55.00   Max.   :81.00

Gambaran awal (Interaksi antara variabel seed dan extract)

# Pastikan seed dan extract sebagai faktor
seeds$seed <- factor(seeds$seed, labels = c("O73", "O75"))
seeds$extract <- factor(seeds$extract, labels = c("Bean", "cucumber"))

seeds

   plate seed  extract  r  n
1      1  O75     Bean 10 39
2      2  O75     Bean 23 62
3      3  O75     Bean 23 81
4      4  O75     Bean 26 51
5      5  O75     Bean 17 39
6      6  O73     Bean  8 16
7      7  O73     Bean 10 30
8      8  O73     Bean  8 28
9      9  O73     Bean 23 45
10    10  O73     Bean  0  4
11    11  O75 cucumber  5  6
12    12  O75 cucumber 53 74
13    13  O75 cucumber 55 72
14    14  O75 cucumber 32 51
15    15  O75 cucumber 46 79
16    16  O75 cucumber 10 13
17    17  O73 cucumber  3 12
18    18  O73 cucumber 22 41
19    19  O73 cucumber 15 30
20    20  O73 cucumber 32 51
21    21  O73 cucumber  3  7

# Buat interaction plot
interaction.plot(x.factor = seeds$seed,  # Faktor pada sumbu X
                 trace.factor = seeds$extract,  # Faktor yang ditampilkan sebagai garis
                 response = seeds$r,  # Variabel respon
                 type = "b",  # Menampilkan titik dan garis
                 pch = c(1, 16),  # Simbol titik (terbuka dan tertutup)
                 col = c("black", "black"),  # Warna garis
                 lty = c(2, 1),  # Jenis garis (putus-putus dan penuh)
                 xlab = "X1", 
                 ylab = "log(r/N)", 
                 legend = TRUE)

Berdasarkan Gambar diatas, dapat diberikan insight bahwa Ekstrak akar mentimun secara signifikan meningkatkan perkecambahan benih aegyptiao 75, sedangkan ekstrak akar kacang tidak memberikan efek yang sama. Analisis regresi logistik (dengan adanya interaksi) bisa digunakan untuk menguji hubungan ini secara statistik.

Selanjutnya, Relabelling untuk memudahkan interpretasi

library(dplyr)

# Relabelling seed dan extract
seeds <- seeds %>%
  mutate(
    seed = ifelse(seed == "O75", 0, 1),        # 0 jika O75, 1 jika O73
    extract = ifelse(extract == "Bean", 0, 1)  # 0 jika Bean, 1 jika Cucumber
  )

# Menampilkan hasil setelah perubahan
seeds

   plate seed extract  r  n
1      1    0       0 10 39
2      2    0       0 23 62
3      3    0       0 23 81
4      4    0       0 26 51
5      5    0       0 17 39
6      6    1       0  8 16
7      7    1       0 10 30
8      8    1       0  8 28
9      9    1       0 23 45
10    10    1       0  0  4
11    11    0       1  5  6
12    12    0       1 53 74
13    13    0       1 55 72
14    14    0       1 32 51
15    15    0       1 46 79
16    16    0       1 10 13
17    17    1       1  3 12
18    18    1       1 22 41
19    19    1       1 15 30
20    20    1       1 32 51
21    21    1       1  3  7

3.2 Model regresi logistik (GLM)

# Model regresi logistik
m1 <- glm(cbind(r, n - r) ~ seed*extract, family = binomial, data = seeds)
summary(m1)


Call:
glm(formula = cbind(r, n - r) ~ seed * extract, family = binomial, 
    data = seeds)

Coefficients:
             Estimate Std. Error z value Pr(>|z|)    
(Intercept)   -0.5582     0.1260  -4.429 9.46e-06 ***
seed           0.1459     0.2232   0.654   0.5132    
extract        1.3182     0.1775   7.428 1.10e-13 ***
seed:extract  -0.7781     0.3064  -2.539   0.0111 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 98.719  on 20  degrees of freedom
Residual deviance: 33.278  on 17  degrees of freedom
AIC: 117.87

Number of Fisher Scoring iterations: 4

Keterangan

extract * seed: Term interaksi antara seed dan extract, sehingga model akan mempertimbangkan efek gabungan keduanya.

Digunakan term `extract * seed:, dikarenakan beberapa pertimbangan, diantaranya:

Dengan extract * seed, model juga memperhitungkan apakah efek extract terhadap perkecambahan berbeda tergantung pada seed yang digunakan/ memastikan bahwa model tidak hanya melihat efek masing-masing variabel secara terpisah, tetapi juga mempertimbangkan bagaimana keduanya berinteraksi dalam mempengaruhi perkecambahan.
Jika hanya menggunakan extract + seed, model hanya memperhitungkan efek masing-masing variabel secara terpisah (efek utama).
Dalam eksperimen ini, ada bukti bahwa cucumber bekerja baik pada satu jenis seed tetapi tidak yang lain, sedangkan bean memiliki efek yang lebih kecil atau berlawanan

Dari output diatas, maka model dapat dituliskan sebagai berikut

Persamaan model

\[ \log \left( \frac{p}{1 - p} \right) = -0.5582 + 0.1459 \cdot \text{seed} + 1.3182 \cdot \text{extract} - 0.7781 \cdot (\text{seed} \times \text{extract}) \]

dimana:

\(\beta_0\) = -0.5582, Intercept, nilai log-odds saat seed = O75 dan extract = Bean.
\(\beta_1 = 0.1459\), Efek perubahan O75 ke O73 pada log-odds.
\(\beta_2 = 1.3182\), Efek perubahan Bean ke Cucumber pada log-odds.
\(\beta_3 = -0.7781\), Efek interaksi antara Seed (O73) dan Extract (Cucumber).

Interpretasi koefisien

1.Intercept (-0.5582)

Saat Seed = O75 dan Extract = Bean, maka log-odds keberhasilan adalah -0.5582.
Jika dikonversi ke probabilitas

\[p=\frac{e^{-0.5582}}{1+e^{-0.5582}} \approx 0.36\]

Jadi probabilitas keberhasilan awalnya sekitar 36%

2.Seed (0.1459, tidak signifikan (\(p = 0.5132\))

Jika Seed berubah dari O75 ke O73, log-odds hanya meningkat 0.1459, tetapi tidak signifikan.
Artinya, jenis benih tidak terlalu berpengaruh secara langsung terhadap probabilitas keberhasilan.

3.Extract (1.3182)

Jika Extract berubah dari Bean ke Cucumber, log-odds meningkat signifikan sebesar 1.3182.
Ini berarti Cucumber meningkatkan peluang keberhasilan secara signifikan dibandingkan Bean.
Dalam probabilitas

\[e^{1.3182} \approx 3.74\]

Jadi, penggunaan Cucumber meningkatkan odds keberhasilan sekitar 3.74 kali lipat dibandingkan Bean.

4.Interaksi Seed dan Extract (\(-0.7781\))

Efek kombinasi Seed = O73 dan Extract = Cucumber mengurangi log-odds sebesar \(0.7781\).
Artinya, meskipun Cucumber secara umum meningkatkan probabilitas keberhasilan, efek ini berkurang jika Seed = O73.

Pembaca ingin sampai prediksi?? Silahkan dicoba sendiri ya.

3.3 Perluasan Model regresi logistik (GLMM)

Selanjutnya, misalkan kita ingin mengetahui apakah terdapat perbedaan sistematis antar plate yang mempengaruhi tingkat keberhasilan germinasi. Dalam hal ini, plate kita anggap sebagai efek acak.

Sehingga kita dapat menggunakan GLMM.

library(lme4)
ns<-21
seeds$s<-factor(1:ns)
m2 <- glmer(cbind(r, n - r) ~ seed*extract + (1 | s), family = binomial(link="logit"), data = seeds)
summary(m2)

Generalized linear mixed model fit by maximum likelihood (Laplace
  Approximation) [glmerMod]
 Family: binomial  ( logit )
Formula: cbind(r, n - r) ~ seed * extract + (1 | s)
   Data: seeds

     AIC      BIC   logLik deviance df.resid 
   117.5    122.8    -53.8    107.5       16 

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-1.60042 -0.78762  0.04326  0.72641  1.24275 

Random effects:
 Groups Name        Variance Std.Dev.
 s      (Intercept) 0.05503  0.2346  
Number of obs: 21, groups:  s, 21

Fixed effects:
             Estimate Std. Error z value Pr(>|z|)    
(Intercept)  -0.54848    0.16608  -3.302 0.000958 ***
seed          0.09743    0.27736   0.351 0.725390    
extract       1.33681    0.23618   5.660 1.51e-08 ***
seed:extract -0.81004    0.38417  -2.109 0.034986 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
            (Intr) seed   extrct
seed        -0.600              
extract     -0.702  0.412       
seed:extrct  0.431 -0.705 -0.617

Berdasarkan hasil GLMM yang Anda berikan, model persamaannya dapat dituliskan sebagai berikut:

\[\text{logit} (p_j) = \log \left( \frac{p_j}{1 - p_j} \right) = \beta_0 + \beta_1 (\text{seed}) + \beta_2 (\text{extract}) + \beta_3 (\text{seed} \times \text{extract}) + u_j\]

atau dalam bentuk nilai estimasi

\[\text{logit} (p_j) = -0.54848 + 0.09743 (\text{seed}) + 1.33681 (\text{extract}) - 0.81004 (\text{seed} \times \text{extract}) + u_j\]

dengan asumsi:

Seed = 0 jika O75, dan Seed=1 jika O73.
extract = 0 jika Bean dan extract = 1 jika Cucumber
\(u_j\) adalah efek acak.

Interpretasi Persamaan

1.Intercept (\(-0.54848\))

Jika tanpa perlakuan (seed = 0 dan extract = 0), maka log-odd rasio keberhasilan adalah \(-0.54848\).
Jika dikonversi ke probailitas.

\[p_j=\frac{e^{-0.54848}}{1+e^{-0.54848}} \approx 0.366\]

Ini berarti ketika tidak ada seed maupun extract, probabilitas keberhasilan sekitar \(36.6\%\).

Efek Seed (\(0.09743\))

Jika seed berubah dari O75(0) -> O73(1), log-odds keberhasilan meningkat sebesar \(0.09743\).
Dalam bentuk probailitas:

\[e^{0.09743} \approx 1.102\]

Artinya, menggunakan O73 dibandingkan O75 meningkatkan odds keberhasilan sebesar \(1.102-1 = 0.102 = 10.2\%\), tetapi \(p-val = 0.725\) menunjukkan bahwa perbedaan ini tidak signifikan secara statistik.

Efek Extract (\(1.33681\))

Jika extract berubah dari Bean(0) -> Cucumber(1), maka log-odds keberhasilan meningkat sebesar \(1.33681\).
Dalam bentuk probabilitas

\[e^{1.33681} \approx 3.81\]

Artinya, menggunakan extract Cucumber meningkatkan odds keberhasilan sekitar 3.81 kali lipat dibandingkan extract Bean, dan efek ini sangat signifikan (\(p < 0.001\)).
Jika hanya ingin meningkatkan keberhasilan, penggunaan extract Cucumber lebih penting dibandingkan perubahan seed dari O75 ke O73.

Interaksi Seed x Extract (\(-0.81004\))

Jika seed = O73 dan extract = Cucumber (keduanya 1), maka terjadi penurunan log-odds keberhasilan sebesar 0.81004 dibandingkan efek masing-masing secara terpisah.
Artinya, kombinasi O73 + Cucumber tidak memberikan keuntungan sebesar yang diperkirakan jika melihat efek seed dan extract secara terpisah.
Dalam bentuk probabilitas

\[e^{-0.81004} \approx 0.44\]
Kombinasi O73 + Cucumber mengurangi odds keberhasilan menjadi 44.4% dari efek masing-masing.
Efek interaksi ini signifikan dengan p-value = \(0.0349\), artinya pengaruhnya cukup penting.

Interpretasi dari output efek acak

Standar deviasi dari efek acak ini adalah 0.2346, yang menunjukkan seberapa besar penyimpangan unit individu dari rata-rata intersept.
Efek acak ini menunjukkan bahwa ada variasi antar-unit yang tidak dapat dijelaskan hanya dengan variabel tetap (fixed effects).
Untuk mengetahui apakah efek acak ini signifikan, bisa dilakukan uji likelihood ratio test (LRT) dengan membandingkan model GLMM dengan model tanpa efek acak, sebagai berikut

library(lme4)

# Model dengan efek acak
m2 <- glmer(cbind(r, n - r) ~ seed*extract + (1 | s), family = binomial(link="logit"), data = seeds)

library(lmtest)
lrtest(m1, m2)

Likelihood ratio test

Model 1: cbind(r, n - r) ~ seed * extract
Model 2: cbind(r, n - r) ~ seed * extract + (1 | s)
  #Df  LogLik Df  Chisq Pr(>Chisq)
1   4 -54.937                     
2   5 -53.770  1 2.3349     0.1265

qchisq(0.95, df = 1) #chi square tabel

[1] 3.841459

Karena \(\chi^2_{hit}= 2.334 < 3.84 = \chi^2_{table}\) dan \(p-val = 0.1265 > 0.05 = \alpha\), maka kita tidak menolak \(H_0\) yang artinya penambahan efek acak (1 | s) dalam model tidak signifikan secara statistik pada tingkat signifikansi \(5\%\).
Dengan kata lain, model dengan efek acak tidak jauh lebih baik daripada model tanpa efek acak.

Efek Acak \(u_j\)

Efek acak menangkap variasi antar kelompok (misalnya, perbedaan antar blok atau individu).
Varians efek acak adalah \(0.05503\), yang berarti ada sedikit variasi antar kelompok tetapi tidak terlalu besar.

Nilai efek acak dapat diekstrak sebagai berikut.

ranef(m2)

$s
    (Intercept)
1  -0.158499495
2   0.009041983
3  -0.182688388
4   0.241332870
5   0.099402378
6   0.080617869
7  -0.066270464
8  -0.117039535
9   0.188888936
10 -0.081427028
11  0.044993839
12  0.062779032
13  0.166013554
14 -0.104336158
15 -0.231904303
16  0.050760258
17 -0.152424836
18  0.025502229
19 -0.022114516
20  0.179552078
21 -0.031746814

with conditional variances for "s"

Nilai dalam tabel diatas adalah penyesuaian (shrinkage) terhadap intercept untuk setiap grup.
Jika efek acak memiliki nilai positif, berarti grup tersebut memiliki intercept yang lebih tinggi dari rata-rata.
Jika efek acak memiliki nilai negatif, berarti grup tersebut memiliki intercept yang lebih rendah dari rata-rata.
Semakin besar (positif atau negatif) nilai efek acak, semakin jauh intercept grup tersebut dari rata-rata populasi.

Misalnya:

Grup 4 memiliki efek acak 0.2413, berarti intercept untuk grup ini lebih tinggi dari rata-rata, sehingga modelnya jika pada plate ke-4 menjadi

\[\text{logit} (p_4) = -0.54848 + 0.09743 (\text{seed}) + 1.33681 (\text{extract}) - 0.81004 (\text{seed} \times \text{extract}) + 0.2413\]

Grup 14 memiliki efek acak -0.1043, berarti intercept untuk grup ini lebih rendah dari rata-rata., sehingga modelnya jika pada plate ke-14 menjadi

\[\text{logit} (p_{14}) = -0.54848 + 0.09743 (\text{seed}) + 1.33681 (\text{extract}) - 0.81004 (\text{seed} \times \text{extract}) -0.1043\]

Efek acak ini mencerminkan bagaimana masing-masing grup menyimpang dari intercept rata-rata model.