Data yang akan digunakan berasal dari Tabel 11.2 dalam Agresti (2002) yang berkaitan dengan penelitian longitudinal. Penelitian bertujuan untuk membandingkan dua jenis obat (“new” vs. “standard”) untuk pengobatan depresi. Variabel respon adalah kondisi depresi (depression) (“normal” vs. “abnormal”). Kategori Y=1 adalah “normal”, sementara Y=0 adalah “abnormal” yang merupakan kategori referensi.
Variabel prediktor yang digunakan adalah:
variabel diagnosis awal, yakni “mild” dan “severe”. Kategori “mild” merupakan kategori referensi.
jenis obat untuk pengobatan depresi, yakni “new” vs. “standard”. Kategori “standard” merupakan kategori referensi.
waktu pemberian obat: 0, 1 dan 2
Berikut adalah tampilan tiga baris pertama.
Dari data, kita melihat bahwa subjek 1 (id = 1) memiliki diagnosis “mild” dan menerima obat “standard” pada waktu 0, 1, dan 2.
URL <- "http://static.lib.virginia.edu/statlab/materials/data/depression.csv"
dat <- read.csv(URL, stringsAsFactors = TRUE)
dat$id <- factor(dat$id)
dat$drug <- relevel(dat$drug, ref = "standard")
head(dat, n = 3)Secara keseluruhan, pasien yang menjadi unit analisis sebanyak 340
length(unique(dat$id))[1] 340Selanjutnya akan dilihat proporsi respons “Normal” (depression == 1) seiring waktu dalam setiap kategori diagnosis dan jenis obat.
Berikut adalah tabel rata-rata berdasarkan diagnosis, obat, dan waktu.
Note: Cara menghitung Proporsi “Normal” yakni misal jika di suatu kelompok ada 10 pasien, dan 7 di antaranya memiliki respons 1 (Normal), maka proporsi “Normal” adalah 7/10 = 0.7.
Data dikelompokkan berdasarkan:
Diagnosis (“mild” atau “severe”)
Obat (“standard” atau “new”)
Waktu (“0”, “1”, “2”)
Output ini mereproduksi Tabel 11.3 dalam Agresti (2002).
Dapat dilihat bahwa proporsi respons “Normal” meningkat seiring waktu untuk keempat kombinasi diagnosis dan pengobatan, dan terjadi peningkatan yang lebih drastis pada obat “new” dibandingkan obat “standard”.
library(magrittr)
with(dat, tapply(depression, list(diagnose, drug, time), mean)) %>%
ftable() %>%
round(2) 0 1 2
mild standard 0.51 0.59 0.68
new 0.53 0.79 0.97
severe standard 0.21 0.28 0.46
new 0.18 0.50 0.83Pemodelan dengan Generalized Estimating Equations (GEE)
Pemodelan GEE digunakan untuk mengestimasi model marginal untuk mengetahui pengaruh diagnosis, obat, dan waktu terhadap variabel respons depresi.
Dalam penelitian ini akan dimasukkan interaksi antara obat dan waktu, untuk menunjukkan bahwa variabel respons bervariasi tergantung pada jenis obat.
Note: Argumen pada GEE dalam `gee()’ yakni:
data → Menentukan data frame yang digunakan.
family → Misalnya binomial untuk data biner.
id → Menentukan faktor pengelompokan. Dalam kasus ini, faktor pengelompokan adalah kolom id dalam data frame, yang menunjukkan setiap individu dalam penelitian.
corstr → Menentukan struktur korelasi antar pengamatan dalam satu individu.
# install.packages("gee")
library(gee)
dep_gee <- gee(depression ~ diagnose + drug*time,
data = dat,
id = id,
family = binomial,
corstr = "independence", scale.fix = TRUE) (Intercept) diagnosesevere drugnew time drugnew:time
-0.02798843 -1.31391092 -0.05960381 0.48241209 1.01744498 summary(dep_gee)
GEE: GENERALIZED LINEAR MODELS FOR DEPENDENT DATA
gee S-function, version 4.13 modified 98/01/27 (1998)
Model:
Link: Logit
Variance to Mean Relation: Binomial
Correlation Structure: Independent
Call:
gee(formula = depression ~ diagnose + drug * time, id = id, data = dat,
family = binomial, corstr = "independence", scale.fix = TRUE)
Summary of Residuals:
Min 1Q Median 3Q Max
-0.94844242 -0.40683252 0.05155758 0.38830952 0.80242231
Coefficients:
Estimate Naive S.E. Naive z Robust S.E. Robust z
(Intercept) -0.02798843 0.1639083 -0.1707566 0.1741865 -0.1606808
diagnosesevere -1.31391092 0.1464151 -8.9738733 0.1459845 -9.0003423
drugnew -0.05960381 0.2222080 -0.2682343 0.2285385 -0.2608042
time 0.48241209 0.1147626 4.2035644 0.1199350 4.0222784
drugnew:time 1.01744498 0.1887954 5.3891398 0.1876938 5.4207709
Estimated Scale Parameter: 1
Number of Iterations: 1
Working Correlation
[,1] [,2] [,3]
[1,] 1 0 0
[2,] 0 1 0
[3,] 0 0 1Koefisien yang ditampilkan berada dalam skala logit, sehingga interpretasi koefisien dilakukan dengan cara yang sama seperti model regresi logistik pada umumnya. * Interpretasi Koefisien Waktu time Koefisien waktu time adalah 0.48. Interpretasi: Kecenderungan pasien yang mengkonsumsi obat standar meningkat sebesar (exp(0.48)=1.61) kali untuk mendapatkan respon “Normal” pada setiap kenaikan satu unit waktu. *Interpretasi Interaksi (drugnew:time) Koefisien interaksi antara obat baru dan waktu (drugnew:time) adalah 1.02. Untuk menghitung efek waktu pada obat baru, kita menambahkan koefisien waktu (0.48) dengan koefisien interaksi (1.02), yakni 0.48 + 1.02 = 1.5
Jika dieksponensialkan (exp(1.5)=4.5). Artinya, kecenderungan pasien yang mengkonsumsi obat baru untuk mendapatkan respons “Normal” 4.5 kali lebih besar untuk setiap kenaikan satu unit waktu.
Pada output tersebut working correlation structure yang digunakan adalah independence, artinya diasumsikan tidak terdapat korelasi antar pengamatan di dalam subjek. Matriks berukuran 3 x 3 karena terdapat 3 observasi per subjek.
Selanjutnya akan dimodelkan jika working correlation structure yang digunakan adalah exchangeable
dep_gee2 <- gee(depression ~ diagnose + drug*time,
data = dat,
id = id,
family = binomial,
corstr = "exchangeable", scale.fix = TRUE) (Intercept) diagnosesevere drugnew time drugnew:time
-0.02798843 -1.31391092 -0.05960381 0.48241209 1.01744498 summary(dep_gee2)
GEE: GENERALIZED LINEAR MODELS FOR DEPENDENT DATA
gee S-function, version 4.13 modified 98/01/27 (1998)
Model:
Link: Logit
Variance to Mean Relation: Binomial
Correlation Structure: Exchangeable
Call:
gee(formula = depression ~ diagnose + drug * time, id = id, data = dat,
family = binomial, corstr = "exchangeable", scale.fix = TRUE)
Summary of Residuals:
Min 1Q Median 3Q Max
-0.94843397 -0.40683122 0.05156603 0.38832332 0.80238627
Coefficients:
Estimate Naive S.E. Naive z Robust S.E. Robust z
(Intercept) -0.02809866 0.1637503 -0.1715945 0.1741791 -0.1613205
diagnosesevere -1.31391033 0.1459325 -9.0035505 0.1459630 -9.0016667
drugnew -0.05926689 0.2221626 -0.2667725 0.2285569 -0.2593091
time 0.48246420 0.1149581 4.1968686 0.1199383 4.0226037
drugnew:time 1.01719312 0.1890913 5.3793750 0.1877014 5.4192084
Estimated Scale Parameter: 1
Number of Iterations: 2
Working Correlation
[,1] [,2] [,3]
[1,] 1.000000000 -0.003432732 -0.003432732
[2,] -0.003432732 1.000000000 -0.003432732
[3,] -0.003432732 -0.003432732 1.000000000Matriks “Working Correlation” menunjukkan estimasi sebesar -0.003433 sebagai korelasi umum antara pasangan observasi dalam satu subjek. Nilai ini relatif sangat kecil, sehingga sebetulnya data ini bisa saja diperlakukan sebagai 1020 observasi independen, bukan 340 subjek dengan masing-masing 3 observasi. Selain itu, estimasi koefisien hampir identik dengan model sebelumnya yang menggunakan struktur korelasi independen.
Struktur korelasi lainnya yang dapat digunakan adalah struktur autoregresif. Struktur ini memungkinkan korelasi antara pengukuran yang diambil dalam waktu yang lebih dekat menjadi lebih tinggi dibandingkan dengan pengukuran yang diambil dalam waktu yang lebih jauh.
Struktur autoregresif yang paling umum adalah AR-1. Dalam penelitian ini berarti bahwa pengukuran pada waktu 1 dan 2 memiliki korelasi tertentu, sedangkan pengukuran pada waktu 1 dan 3 memiliki korelasi yang lebih kecil.
dep_gee3 <- gee(depression ~ diagnose + drug*time,
data = dat,
id = id,
family = binomial,
corstr = "AR-M", Mv = 1, scale.fix = TRUE) (Intercept) diagnosesevere drugnew time drugnew:time
-0.02798843 -1.31391092 -0.05960381 0.48241209 1.01744498 summary(dep_gee3)
GEE: GENERALIZED LINEAR MODELS FOR DEPENDENT DATA
gee S-function, version 4.13 modified 98/01/27 (1998)
Model:
Link: Logit
Variance to Mean Relation: Binomial
Correlation Structure: AR-M , M = 1
Call:
gee(formula = depression ~ diagnose + drug * time, id = id, data = dat,
family = binomial, corstr = "AR-M", Mv = 1, scale.fix = TRUE)
Summary of Residuals:
Min 1Q Median 3Q Max
-0.94844464 -0.40691023 0.05155536 0.38824284 0.80236892
Coefficients:
Estimate Naive S.E. Naive z Robust S.E. Robust z
(Intercept) -0.02770314 0.1643892 -0.1685216 0.1741163 -0.1591071
diagnosesevere -1.31386512 0.1472418 -8.9231818 0.1459894 -8.9997314
drugnew -0.05959816 0.2226449 -0.2676826 0.2284592 -0.2608700
time 0.48240753 0.1147578 4.2036998 0.1199244 4.0225960
drugnew:time 1.01732678 0.1887941 5.3885527 0.1876727 5.4207510
Estimated Scale Parameter: 1
Number of Iterations: 2
Working Correlation
[,1] [,2] [,3]
[1,] 1.000000e+00 0.008477443 7.186704e-05
[2,] 8.477443e-03 1.000000000 8.477443e-03
[3,] 7.186704e-05 0.008477443 1.000000e+00Dari hasil tersebut “Working Correlation” menunjukkan nilai yang sangat kecil. Parameter korelasi diestimasi sebesar 0.008477 (korelasi antara pengukuran 1-2 dan 2-3).
Salah satu struktur korelasi lain adalah matriks unstructured. Struktur ini memungkinkan semua korelasi bervariasi secara bebas.
dep_gee4 <- gee(depression ~ diagnose + drug*time,
data = dat,
id = id,
family = binomial,
corstr = "unstructured", scale.fix = TRUE) (Intercept) diagnosesevere drugnew time drugnew:time
-0.02798843 -1.31391092 -0.05960381 0.48241209 1.01744498 summary(dep_gee4)
GEE: GENERALIZED LINEAR MODELS FOR DEPENDENT DATA
gee S-function, version 4.13 modified 98/01/27 (1998)
Model:
Link: Logit
Variance to Mean Relation: Binomial
Correlation Structure: Unstructured
Call:
gee(formula = depression ~ diagnose + drug * time, id = id, data = dat,
family = binomial, corstr = "unstructured", scale.fix = TRUE)
Summary of Residuals:
Min 1Q Median 3Q Max
-0.94773674 -0.40645713 0.05226326 0.38927858 0.79975454
Coefficients:
Estimate Naive S.E. Naive z Robust S.E. Robust z
(Intercept) -0.02552611 0.1679741 -0.1519645 0.1726392 -0.1478581
diagnosesevere -1.30484850 0.1461691 -8.9269772 0.1450136 -8.9981088
drugnew -0.05438636 0.2282121 -0.2383149 0.2271321 -0.2394481
time 0.47587182 0.1160832 4.0994035 0.1190418 3.9975178
drugnew:time 1.01297603 0.1887379 5.3671034 0.1865407 5.4303205
Estimated Scale Parameter: 1
Number of Iterations: 3
Working Correlation
[,1] [,2] [,3]
[1,] 1.00000000 0.07393977 -0.02741128
[2,] 0.07393977 1.00000000 -0.05669559
[3,] -0.02741128 -0.05669559 1.00000000Bagaimana cara memilih struktur korelasi yang tepat? Menurut Agresti (2002), estimasi GEE tetap valid meskipun struktur korelasi yang digunakan tidak tepat. Agresti menyarankan untuk memulai dengan struktur exchangeable, lalu memeriksa bagaimana estimasi koefisien dan standar error berubah ketika menggunakan struktur korelasi lain. Jika perubahannya minimal, maka pilihlah struktur korelasi yang lebih sederhana.
Sumber:
Agresti, A. (2002). Categorical data analysis (2nd ed.). Wiley
https://library.virginia.edu/data/articles/getting-started-with-generalized-estimating-equations
#Sebagai contoh kasus, simulasi pemodelan untuk memprediksi jumlah kunjungan wisatawan di sebuah kawasan wisata alam (seperti kawasan Danau Toba).Dalam analisis pariwisata, faktor-faktor prediktor sering kali saling berkorelasi sangat tinggi (multikolinearitas). Misalnya, daerah yang memiliki "Jumlah Hotel" banyak hampir pasti juga memiliki "Jumlah Restoran" yang banyak dan "Anggaran Promosi" yang tinggi. Jika menggunakan regresi linier biasa (MLE), multikolinearitas ini akan membuat estimasi koefisien menjadi tidak stabil. Di sinilah Regresi Ridge (dengan penalti L2) menyelamatkan model.
# 1. Memuat library yang dibutuhkan
library(glmnet)
# 2. Simulasi Data Pariwisata
set.seed(123)
n <- 100 # Jumlah observasi (misal: data bulanan atau titik lokasi pengamatan)
# Membuat prediktor yang saling berkorelasi tinggi (masalah multikolinearitas)
jumlah_hotel <- rnorm(n, mean = 50, sd = 10)
# Restoran dan promosi berkorelasi linier dengan jumlah hotel
jumlah_restoran <- jumlah_hotel * 1.5 + rnorm(n, mean = 0, sd = 5)
anggaran_promosi <- jumlah_hotel * 0.8 + rnorm(n, mean = 0, sd = 10)
# Variabel independen lainnya
kualitas_infrastruktur <- rnorm(n, mean = 7, sd = 1)
indeks_cuaca <- rnorm(n, mean = 20, sd = 5)
# Variabel Respon: Jumlah Kunjungan Wisatawan
kunjungan_wisatawan <- 100 + 5*jumlah_hotel + 3*jumlah_restoran +
2*anggaran_promosi + 10*kualitas_infrastruktur +
0.5*indeks_cuaca + rnorm(n, mean = 0, sd = 50)
# Menyatukan ke dalam Data Frame
data_wisata <- data.frame(kunjungan_wisatawan, jumlah_hotel, jumlah_restoran,
anggaran_promosi, kualitas_infrastruktur, indeks_cuaca)
# 3. Persiapan Data untuk glmnet
# glmnet tidak menerima data frame secara langsung.
# Kita harus memisahkannya menjadi matriks prediktor (x) dan vektor respon (y).
x <- model.matrix(kunjungan_wisatawan ~ ., data = data_wisata)[, -1] # [, -1] membuang kolom intercept bawaan
y <- data_wisata$kunjungan_wisatawan
# 4. Membangun Model Regresi Ridge
# Argumen alpha = 0 adalah kunci untuk memberitahu glmnet agar menggunakan penalti L2 (Ridge)
model_ridge <- glmnet(x, y, alpha = 0)
# 5. Mencari Nilai Lambda Optimal dengan Cross-Validation
# Kita membiarkan data menentukan seberapa besar penalti yang paling pas
cv_ridge <- cv.glmnet(x, y, alpha = 0)
lambda_optimal <- cv_ridge$lambda.min
cat("Nilai Lambda Optimal yang ditemukan:", lambda_optimal, "\n\n")Nilai Lambda Optimal yang ditemukan: 9.537566 # 6. Menampilkan Koefisien Akhir
# Menggunakan lambda optimal untuk melihat koefisien yang sudah disusutkan (shrinkage)
koefisien_ridge <- predict(model_ridge, s = lambda_optimal, type = "coefficients")
print("Koefisien Regresi Ridge:")[1] "Koefisien Regresi Ridge:"print(koefisien_ridge)6 x 1 sparse Matrix of class "dgCMatrix"
s=9.537566
(Intercept) 120.004763
jumlah_hotel 3.827045
jumlah_restoran 3.462542
anggaran_promosi 2.023730
kualitas_infrastruktur 5.355834
indeks_cuaca 2.149608# 1. Memuat Library
library(glmnet)
# 2. Simulasi Data Transportasi Publik
set.seed(42)
n <- 200 # Jumlah area pengamatan
# Membuat variabel tata kota yang saling berkorelasi (Multikolinearitas)
kepadatan_penduduk <- rnorm(n, mean = 8000, sd = 2000)
jumlah_halte <- kepadatan_penduduk * 0.005 + rnorm(n, mean = 0, sd = 10)
jarak_pusat_kota <- kepadatan_penduduk * -0.002 + rnorm(n, mean = 20, sd = 5)
# Variabel Respon: Jumlah penumpang bus
penumpang <- 100 + 0.05*kepadatan_penduduk + 15*jumlah_halte - 8*jarak_pusat_kota + rnorm(n, mean = 0, sd = 300)
# Menyiapkan matriks X dan vektor y
X <- cbind(kepadatan_penduduk, jumlah_halte, jarak_pusat_kota)
y <- penumpang
# 3. Fitting Model Regresi Ridge
# Mencari lambda optimal menggunakan cross-validation
cv_ridge <- cv.glmnet(X, y, alpha = 0)
best_lambda <- cv_ridge$lambda.min
# Membangun model final dengan lambda terbaik
model_ridge <- glmnet(X, y, alpha = 0, lambda = best_lambda)
# 4. MENGHITUNG LOG-LIKELIHOOD MANUALLY
# A. Dapatkan prediksi (y-hat) dan residual (error)
y_hat <- predict(model_ridge, newx = X)
residual <- y - y_hat
# B. Hitung Residual Sum of Squares (RSS) dan Varians (sigma kuadrat)
rss <- sum(residual^2)
sigma2 <- rss / n # Menggunakan pembagi n untuk estimator Maximum Likelihood
# C. Hitung Log-Likelihood Klasik (Tanpa Penalti)
log_lik_klasik <- -0.5 * n * log(2 * pi) - 0.5 * n * log(sigma2) - (1 / (2 * sigma2)) * rss
# D. Hitung Suku Penalti L2 (Ridge)
# Ekstrak koefisien (kecuali intercept)
beta_coef <- coef(model_ridge)[-1, ]
penalti_L2 <- best_lambda * sum(beta_coef^2)
# E. Hitung Penalized Log-Likelihood
penalized_log_lik <- log_lik_klasik - penalti_L2
# 5. Menampilkan Hasil
cat("--- HASIL KALKULASI LOG-LIKELIHOOD ---\n")--- HASIL KALKULASI LOG-LIKELIHOOD ---cat("Nilai Lambda Optimal :", best_lambda, "\n")Nilai Lambda Optimal : 27.16864 cat("Log-Likelihood Klasik :", log_lik_klasik, "\n")Log-Likelihood Klasik : -1409.16 cat("Besaran Penalti L2 :", penalti_L2, "\n")Besaran Penalti L2 : 7060.706 cat("Penalized Log-Likelihood:", penalized_log_lik, "\n")Penalized Log-Likelihood: -8469.866